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Electronic excitations lie at the origin of most of the commonly measured spectra. However, the 
first-principles computation of excited states requires a larger effort than ground-state calculations, 
which can be very efficiently carried out within density-functional theory. On the other hand, two 
theoretical and computational tools have come to prominence for the description of electronic 
excitations. One of them, many-body perturbation theory, is based on a set of Green's-function 
equations, starting with a one-electron propagator and considering the electron-hole Green's function 
for the response. Key ingredients are the electron's self-energy 2 and the electron-hole interaction. A 
good approximation for 2 is obtained with Hedin's GW approach, using density-functional theory as 
a zero-order solution. First-principles GW calculations for real systems have been successfully carried 
out since the 1980s. Similarly, the electron-hole interaction is well described by the Bethe-Salpeter 
equation, via a functional derivative of 2. An alternative approach to calculating electronic excitations 
is the time-dependent density-functional theory (TDDFT), which offers the important practical 
advantage of a dependence on density rather than on multivariable Green's functions. This approach 
leads to a screening equation similar to the Bethe-Salpeter one, but with a two-point, rather than a 
four-point, interaction kernel. At present, the simple adiabatic local-density approximation has given 
promising results for finite systems, but has significant deficiencies in the description of absorption 
spectra in solids, leading to wrong excitation energies, the absence of bound excitonic states, and 
appreciable distortions of the spectral line shapes. The search for improved TDDFT potentials and 
kernels is hence a subject of increasing interest. It can be addressed within the framework of 
many-body perturbation theory: in fact, both the Green's functions and the TDDFT approaches profit 
from mutual insight. This review compares the theoretical and practical aspects of the two approaches 
and their specific numerical implementations, and presents an overview of accomplishments and work 
in progress. 
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I. INTRODUCTION 

In every spectroscopic experiment one perturbs the 
sample (by incoming photons, electrons, etc.) and mea- 
sures the response of the system to this perturbation. In 
other words, the system is excited. Therefore it is in gen- 
eral not sufficient to calculate ground-state properties in 
order to interpret or predict results of experiments like 
photoemission, electron-energy loss, absorption, etc. 

Direct and inverse photoemission and absorption can 
be taken as the prototype spectroscopies which one 
would like to describe in this context. They are sche- 
matically depicted in Fig. 1. In the photoemission pro- 
cess the system absorbs a photon h v, and an electron is 
ejected whose kinetic energy E k is then measured at 
some distance. If one considers this photoelectron to be 
completely decoupled from the sample, energy and mo- 
mentum conservation allow one to deduce the change in 
total energy of the sample, which is interpreted as the 
energy level of the "hole," i.e., the level that was for- 
merly occupied by the photoelectron. Hence as a first 
approach one can state that photoemission measures the 
density of occupied states. By analogy, inverse photo- 
emission yields information about the density of unoc- 
cupied states. In absorption experiments, an electron is 
excited from an occupied state into a conduction state. 
This process looks, at first glance, like the sum of a pho- 



toemission and an inverse photoemission experiment 
(creation of a hole and an electron). Instead, as shown in 
Fig. 1, in the absorption measurement the excited elec- 
tron remains inside the system and cannot be supposed 
to be a free electron decoupled from the others. Hence 
whereas direct and inverse photoemission results are of- 
ten already well described by the density of occupied 
and unoccupied states, respectively, one realizes that (i) 
in absorption, the joint density of occupied and unoccu- 
pied states must be considered, (ii) even over a small 
range of excitation energies, transition probabilities can 
vary considerably and must be taken into account, and 
(hi) the excited electron and the hole cannot be treated 
separately, since the electron feels the presence of the 
hole. Point (iii) constitutes the main difficulty for a cor- 
rect description of this type of experiment. 

The interpretation of photoemission spectra as a den- 
sity of occupied states is linked to the picture of inde- 
pendent electrons which occupy some well-defined en- 
ergy level in the system. Of course, electrons are not 
independent, and it is clear that, for example, an elec- 
tron that leaves the sample will lead the remaining elec- 
trons to relax. This relaxation energy and other 
quantum-mechanical contributions must be taken into 
account if energy differences are to be calculated cor- 
rectly. In other words, in photoemission the single- 
electron energy levels are renormalized by the presence 
of the other electrons. One can then still retain the pic- 
ture of one-particle energy levels, but these particles are 
quasielectrons and quasiholes, i.e., they contain the ef- 
fects of all the other particles (Landau, 1957a, 1957b, 
1959). They can still be described by a sort of one- 
particle Schrodinger equation, which does, however, 
contain rather complicated effective potentials (also 
called optical potentials) reflecting these interactions. 
Moreover, it is clear that, in the case of absorption, even 
a very sophisticated one-quasiparticle Schrodinger equa- 
tion would be inadequate, since the quasielectron and 
quasihole must be described simultaneously, requiring 



Direct Inverse Absorption 

photoemission photoemission 




FIG. 1. Schematic representation of the excitations involved in 
direct photoemission, inverse photoemission, and absorption 
spectroscopies. Photoemission can be resolved in angle, spin, 
and time; absorption can be resolved in polarization and time. 
This allows a direct probing of the electronic and structural 
properties of bulk and low-dimensional samples including dy- 
namical effects. \E = Ef— E'f , apart from phonons and radia- 
tive losses. 
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an effective two-particle equation. Progress in the theo- 
retical description of spectroscopy is in fact very often 
linked to progress in finding better effective one- or two- 
particle Hamiltonians. 

One important motivation for seeking more precise 
theoretical descriptions is, of course, the fact that dis- 
crepancies arise between experimental spectra and the 
spectra calculated at some approximate level. This mo- 
tivation becomes stronger as more precise experiments 
become available. 

A. Motivation from experiments 

The excitations considered in the present work con- 
cern electrons, mainly in the valence energy region. 
Most of the theoretical tools presented in the following 
are general enough to be used also for the study of ex- 
citations involving core electrons, although in that case 
one is frequently interested in more complex phenom- 
ena such as the Auger effect (Verdozzi et al, 2001), 
which is not explicitly treated in this review. Traditional 
techniques for the study of valence electron excitations 
use either photons (visible and ultraviolet absorption, 
transmission and reflectivity spectra), electrons 
(electron-energy-loss spectra), or both (electron photo- 
emission and inverse photoemission spectra). 

During the last two decades, important improvements 
in many of these techniques have been achieved, due to 
the ongoing substantial progress in obtaining (a) high 
spectral and spatial resolution, brightness; (b) short 
measurement times (scale of femtoseconds); (c) high 
spatial coherence; and (d) low temperatures (Smith, 
2001). Large contributions come from the progress 
made at synchrotron-radiation sources and with ultrafast 
lasers. In particular, photon beams in the soft-x-ray en- 
ergy region (extreme ultraviolet) are increasingly used, 
due to the availability of high-brilliance sources (Smith, 
2001). Several striking examples can be found in the re- 
cent literature. Femtosecond lasers yield information 
about elementary electronic processes occurring at sur- 
faces on time scales from pico- to femtoseconds, which 
are relevant for potential technological applications. For 
example, an electronic excitation is the initial step in 
many chemical reactions, and the energetics and lifetime 
of this process directly govern the reaction probability. 
Hence chemical selectivity can be obtained through an 
activation of the desired reaction via a femtosecond 
electronic excitation (Sundstrom, 1996; see Diau etal, 
1998, for an example). Time -resolved femtosecond pho- 
toemission spectroscopy has been used to gain insight 
into electronically induced adsorbate reactions at sur- 
faces and kinetics of growth, and to monitor in real time 
and with atomic resolution the dynamics of electrons 
(Plummer, 1997) and the motion of atoms at surfaces 
(Petek etal., 2000). Recent advances in ultrafast elec- 
tron diffraction have yielded direct imaging of transient 
structures in chemical reactions (Ihee et al, 2001). X-ray 
microscopy allowed the study of biological matter, like 
the internal structure of a cell, with a spatial resolution 
of 36 nm using photons in the "water window" energy 



range, i.e., between 290 and 530 eV, where carbon ab- 
sorbs (excitation of the Is electrons) and water is rela- 
tively transparent. Photoelectron spectra with high spa- 
tial resolution are at the basis of photoelectron 
microscopy, where a 20-nm resolution can be obtained 
(Nolting etal., 2000). The determination of the Fermi 
surface of very complex materials such as cuprate super- 
conductors (Zhou et al, 1999), the study of a Fermi gap 
opening in complex structures (Cepek et al, 2001), and 
the measurement of the energy dependence of the line- 
width associated with photoemission peaks near the 
Fermi energy (Valla etal, 1999a, 1999b) are now pos- 
sible because of the improved spectral resolution. 

Often, it is clear that a simple one-particle picture is 
intrinsically inadequate for describing the processes oc- 
curring in the experiments. A good example is resonant 
photoemission, where core-valence absorption and va- 
lence electron Auger emission interfere. Many aspects 
of these experiments can in principle be described by the 
theoretical tools discussed in this review. In the follow- 
ing, we define the common framework of the ap- 
proaches and then consider and compare them in detail. 

B. Theoretical framework 

1. Overview 

This review treats approaches using or leading to the 
picture of effective particles, i.e., "quasiparticles," as 
outlined above. There are, of course, other ways to treat 
the many-body problem, e.g., the configuration- 
interaction approach of quantum chemistry (see, for ex- 
ample, Jensen, 1999; Szabo and Ostlund, 1983; Bonacic- 
Koutecky etal, 1990 for an application). The latter is 
based on the minimization of the energy with respect to 
the expansion coefficients of a trial many-body wave 
function, written as a linear combination of determi- 
nants. A description of configuration-interaction-based 
techniques is clearly beyond the scope of this paper, al- 
though they can be very efficient (e.g., in the determina- 
tion of higher excited states for point defects in Si0 2 ; 
Raghavachari etal, 2002). However, the drawback of 
the configuration-interaction method with respect to 
density-functional-based or Green's-function-based 
"quasiparticle" methods is its unfavorable scaling with 
the system's size. 1 



In particular, full configuration-interaction calculations in- 
volve a number of determinants which grows exponentially 
with the size of the system. Simplified methods have been de- 
vised, such as the truncated configuration interaction, which 
reduces the scaling to M 6 (where M is the number of basis 
functions) when only determinants with a number of excited 
electrons lower than or equal to 2 are included. Other scalings, 
namely, M 8 or M w , are obtained by also including triply or 
quadruply excited determinants. The approaches which will be 
discussed in this review can exhibit a much better scaling, ac- 
cording to the system and to the numerical implementation 
(Benedict et al, 1998a; Bertsch et al, 2000; Hahn et al, 2002; 
Vasiliev et al, 2002). 
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The concept of quasiparticles is intimately linked to 
that of band structure (or discrete energy levels in a 
finite system), and the corresponding one-electron-like 
Schrodinger equations. Band-structure equations can 
contain the complications of electron-electron interac- 
tions implicitly in empirical or semiempirical Hamilto- 
nians, as in the tight-binding approach, or they can treat 
those effects explicitly. It is this latter point which is of 
interest here. We shall hence discuss the quasiparticle 
formalism of many-body perturbation theory. In addi- 
tion, we shall also discuss approaches based on the static 
and time-dependent density-functional theory (DFT and 
TDDFT). 

Static DFT (Hohenberg and Kohn, 1964) in the for- 
mulation of Kohn and Sham (1965) is a ground-state 
theory in the form of an effective one-particle Schro- 
dinger equation (the Kohn-Sham equation). The Kohn- 
Sham eigenvalues are often interpreted as quasiparticle 
energies, and their differences as optical excitation ener- 
gies, without formal justification. The Kohn-Sham eigen- 
values and eigenfunctions are also often used as starting 
point for further excited-state calculations. On the other 
hand, DFT has been extended to time-dependent DFT 
(Runge and Gross, 1984), which is in principle an exact 
theory for the description of neutral excitations (such as 
those involved in absorption). In fact, TDDFT can be 
written in the form of an effective two-particle equation, 
which can be directly compared to the effective two- 
quasiparticle equation of many-body perturbation 
theory. 

In the last few years, several reviews have addressed 
different aspects of the calculation of electronic excita- 
tions in both finite and infinite systems using either the 
quasiparticle picture (see, e.g., Aryasetiawan and Gun- 
narsson, 1998; Strinati 1988; Aulbur etal, 1999; Farid, 
1999a; Hedin, 1999; Rohlfing and Louie, 2000) or TD- 
DFT (see, for example, Gross et al. , 1994, 1996; Casida 
1995, 1996; Dobson, Vignale, and Das, 1997; Rubio 
etal, 1997; van Leeuwen, 2001; Burke etal, 2002). It 
emerges that, although both methods are in principle 
exact when applied to the appropriate problem, they 
present different drawbacks: in the context of many- 
body perturbation theory, the effective one- and two- 
quasiparticle Hamiltonians are essentially well estab- 
lished, but the methods become numerically 
impracticable for complex systems. TDDFT, on the 
other hand, could in principle lead to technically simpler 
equations, but its large-scale application is at present 
prevented by the fact that a good general approximation 
for the corresponding effective two-particle Hamiltonian 
(or equivalent quantities) has not yet been found. The 
large increase of interest in theoretical spectroscopic 
analysis from first principles suggests a joint venture in 
exploring both approaches. The aim is to reach an effec- 
tive method of handling excitations of many-electron 
systems with an effort comparable to that of DFT-based 
approaches used for the calculation of ground-state 
properties. 

In the present work we review some of the essential 
physics contained in the various approximations used to 



describe spectroscopic measurements from first prin- 
ciples. We address in detail some frequently asked ques- 
tions related to the calculation of photoemission and op- 
tical spectra, including those which remain open and 
may be solved in the near future. We attempt to give a 
general perspective of our present understanding of dif- 
ferent spectroscopies and to compare the many-body 
perturbation theory approaches — in particular those 
based on Hedin's equations (Hedin, 1965) — with time- 
dependent density-functional theory. This comparison 
with an eye to practical applications, together with the 
focus on the calculation of photon absorption and 
electron-energy-loss spectra which involve two- 
quasiparticles (electron-hole) excitations, characterizes 
the present work and differentiates it from other re- 
cently published reviews on Green's-function methods, 
more focused on one -particle Green's functions and 
single-quasiparticle excitations, such as those involved in 
photoemission or inverse photoemission (e.g., Aryase- 
tiawan and Gunnarsson, 1998; Aulbur et al, 1999; Farid, 
1999a; Hedin, 1999). 

Since this field of research is rather broad, and in or- 
der to make the discussions more focused and specific, 
we have restricted the review to the description of 
purely electronic excitations from valence states, for 
nonmagnetic systems. We have hence excluded the ex- 
plicit discussion of core-level spectroscopy (Almbladh 
and Hedin, 1983) and Auger spectroscopy (Verdozzi 
et al, 2001), although the theory presented here remains 
valid in that energy range. We also exclude the problem 
of strong correlation and the methods that have been 
developed to treat these problems in particular, e.g., the 
local-density approximation plus an on-site Hubbard re- 
pulsion method (LDA+U; Anisimov etal, 1997) and 
the dynamical mean-field theory (Georges et al, 1996). 

Unless otherwise stated, we shall use atomic units 
throughout the paper (i.e., e 2 = h = m e =l). 

2. Effective Hamiltonians and effective interactions 

The earliest attempts to cast the many-body problem 
into the form of one particle moving in some mean (or, 
more generally, effective) field due to the electron- 
electron interaction are the Hartree (1928) and the 
Hartree-Fock (Fock, 1930) approaches: a product or de- 
terminant ansatz for the many-body wave function al- 
lows one to minimize the total energy variationally, lead- 
ing to an effective one-particle Schrodinger equation. It 
is interesting to note that the Hartree-Fock method al- 
ready yields good results, like total energy differences, 
in certain systems. However, the interpretation of the 
Hartree-Fock eigenvalues as electron addition and re- 
moval energies does not lead to satisfactory agreement 
between theory and experiment. As a guideline we may 
take the minimum direct quasiparticle gap at T in dia- 
mond, for which Mauger and Lannoo (1977) obtained a 
value of E g — 15 eV in a self-consistent Hartree-Fock 
calculation, instead of the experimental quasiparticle 
gap of 7.3 eV (inferred from optical experiments by 
Roberts and Walker, 1967). 
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The need for going beyond Hartree-Fock in electronic 
structure calculations was in fact recognized early. 
Many-body correlation effects have been studied since 
the 1930s (Wigner, 1934); an important milestone is the 
work of Quinn and Ferrell (1958) on the homogeneous 
electron gas and metals. A key role is played by the 
electron self-energy, which is the proper exchange- 
correlation potential acting on an excited electron or 
hole. A sound basis for modern band-structure calcula- 
tions was laid by Hedin (Hedin, 1965; Hedin and Lun- 
dqvist, 1969) with the so-called GW approximation for 
the self-energy, which was applied to the electron gas in 
that work. The essential improvement offered by the 
GW one-particle Schrodinger equation over the 
Hartree-Fock equation lies in the explicit description of 
the potential induced by the additional particle. In order 
to understand this, one might ask why the eigenvalues of 
the Hartree-Fock equation do not yield a satisfactory 
description of quantities like the band structure, al- 
though Koopman's theorem says that it is justified to 
identify eigenvalues with total energy differences. So 
what is missing, and which ingredients should reason- 
ably appear in a theory going beyond Hartree-Fock and 
allowing one to use the eigenvalues of some effective 
Hamiltonian? 

First, Koopman's theorem supposes that the one- 
electron orbitals are frozen upon changing the number 
of electrons. In order to go beyond this drastic approxi- 
mation, one has to include the relaxation of the orbitals 
as a response to the addition of an electron or a hole to 
the system. This effect is contained in the so-called delta- 
self-consistent-field calculations for finite systems, where 
the difference in total energy between two self- 
consistent calculations, for N and jV±1 electrons, is ob- 
tained explicitly. In other words, total energy differences 
are calculated successfully within Hartree-Fock, because 
in the calculation for N± 1 electrons the wave functions 
are allowed to relax with respect to those of the 
jV-electron calculation. This response of the system to 
the additional electron or hole — the screening of the ad- 
ditional particle — can be reformulated in linear response 
in terms of the dielectric matrix e, defined by the relation 

V tot (r)= J dt' e -\r,r')V ext {r') (1.1) 

between the external potential V ext and the screened 
(total) potential V tot . The total potential is the sum of 
the external potential and the potential due to the po- 
larization of the system induced by the external pertur- 
bation. It is the relaxation of the wave functions which 
gives rise to this polarization, and one can therefore try 
to use the concept of the dielectric matrix in order to 
correct the Hartree-Fock one-particle Schrodinger equa- 
tion, without calculating the relaxation explicitly as in 
the delta-self-consistent-field approach. One has hence 
to calculate the induced potential V ind =V tot —V ext due 
to an additional electron, and acting on that electron 
itself. As a first step, one could assume that the addi- 
tional electron in a state n gives rise to an induced po- 
tential that is proportional to the perturbing "charge 



density" | ip n {r)\ 2 (as in the familiar picture of the image 
charge). However, whereas Hartree relaxation effects 
are very important in small systems, in solids in a Bloch 
picture they are negligible, and hence this effect alone 
cannot lead to a satisfactory correction of the Hartree- 
Fock equation. Rather, one should take [^(r)) 2 as the 
probability to find the additional electron in some point 
r, and then, supposing that the additional electron is at r 
(i.e., taking correlation into account), calculate the in- 
duced potential in the same point. Then, in the corre- 
sponding equations | if/ n \ 2 is replaced by a S function, and 
the potential which should be added to the Hartree- 
Fock equation turns out to be the so-called Coulomb- 
hole term \v(e~ l — 1), v being the bare Coulomb 
interaction. 2 

Of course, the change in the charge distribution due 
to the additional, perturbing, electron will also affect the 
exchange term of the original Hartree-Fock equation, 
which now turns out to be screened by e ~ 1 . These argu- 
ments lead to an effective one-particle Schrodinger 
equation for an additional electron (or hole) which 
reads 

+ \ j dr'v(r',r)[ s - l (r, r ')-S(r,r')]j^ n ( r ) 

occ 

- dr'W(r,r')X #(*')&(*) WO = «»&,(*), 

J s=l 

(1.2) 

where W=e~ 1 v is the screened Coulomb interaction. 
This is the equation known as the static COHSEX (Cou- 
lomb hole plus screened exchange) approximation (He- 
din, 1965), where the terms added to the kinetic-energy 
operator and the Hartree potential constitute the self- 
energy ^ COHSEX _ p or diamond, Hybertsen and Louie 
(1986) have found that the energy difference between 
valence and conduction bands is overestimated in the 
COHSEX approximation by about 1 eV, which means, 
however, that most of the error of the Hartree-Fock re- 
sult cited above is removed by this correction. 

To go further, one must take into account the fact that 
the response of a system to an external perturbation is 
frequency dependent, which means that one has to in- 
troduce the concept of a response in time (dynamical 
correlation and memory effects). In fact, the COHSEX 
potential is nothing other than the static limit of the 
exchange-correlation self-energy calculated in the so- 



2 The factor \ describes the fact that a charge e is taken from 
infinity to a point r of the system or, equivalently, built up from 
e = 0. The part of the work coming from the induced charge 
which is needed to do this is J e dqV ind (q) = f e dqq(E~ 1 — l)v 
= l/2e 2 *(W— v). This is nothing other than the adiabatic 
building up of charge density (Hedin, 1965, 1999) of classical 
electrostatics. 
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called GW approximation (Hedin, 1965), which, as 
pointed out above, can be considered as the state-of-the- 
art tool for band-structure calculations today. 

After its introduction for metals, the GW approxima- 
tion was also used early for semiconductors and insula- 
tors, first in static COHSEX (Brinkman and Goodman, 
1966; Lipari and Fowler, 1970; Brener, 1975a, 1975b), 
later using a plasmon-pole approximation in a model 
e _1 (Bennett and Inkson, 1977; Inkson and Bennett, 
1978) and, starting from a Hartree-Fock calculation, us- 
ing a realistic frequency- and wave-vector-dependent di- 
electric matrix, in a basis of localized orbitals (Strinati 
et al, 1980, 1982). 

It was shown that the GW approach managed to cor- 
rect the largest part of the band-gap error; in fact, for 
the example of the gap at T in diamond, Strinati et al. 
(1980) found a value of 7.4 eV, which is half of the 
Hartree-Fock gap and very close to the experimental 
one. 3 The fact that the inclusion of dynamical effects in 
the response functions leads to a reduction of the quasi- 
particle gap with respect to the static COHSEX result is 
a general finding, and the example of diamond shows 
the typical order of magnitude of the dynamical effects 
(i.e., about 10-20 % of the experimental gap). 

On the other hand, G W calculations are cumbersome 
as compared, for example, to DFT ones. As pointed out 
above, the Kohn-Sham equation also has the form of an 
effective one-particle Schrodinger equation, and it is 
therefore tempting to use its eigenvalues as electron ad- 
dition and removal energies. However, this yields large 
errors (Lundqvist and March, 1983); in particular one 
finds a strong underestimate of the band gap of semicon- 
ductors and insulators. For diamond, the difference of 
Kohn-Sham eigenvalues calculated in the local-density 
approximation (LDA) has yielded a gap at T of E g 
= 5.51 eV (Hybertsen and Louie, 1985). The origin of 
this failure, whether it is mainly due to the wrong inter- 
pretation of the Kohn-Sham equation or to the LDA, is 
still discussed today and will be studied in Sec. IV.A.2. In 
any case, DFT is a good starting point for further calcu- 
lations, and starting from the work of Hybertsen and 
Louie (1985, 1986), and Godby et al. (1986, 1987, 1988), 
today ab initio GW calculations are mostly performed 
using Kohn-Sham results as ingredients (to be precise, 
the Kohn-Sham eigenvalues and eigenfunctions are used 
to construct the self-energy of the GW form). These 
fully ab initio GW calculations generally yield very good 
agreement, most often better than 10-15 %, between 
the experimental and the calculated band structure, 
apart from the case of strongly correlated systems. 
Again for the example of diamond, Hybertsen and 
Louie (1985) and Godby et al. (1987) have performed ab 
initio GW calculations and obtained a quasiparticle gap 
at T of £ =7.38 eV and £ ? = 7.26 eV, respectively, 
which compare well with the experimental value of 7.3 
eV (Roberts and Walker, 1967). 



The GW approach thus yields very gratifying results 
concerning the band structure, but having a good band 
structure is not enough when one is interested in spec- 
troscopies like absorption, where one creates a neutral, 
i.e., electron-hole type of excitation. One is in that case 
talking about the response of the system to a (time- 
dependent) external potential, which implies that one is 
interested in the dielectric matrix itself, as defined in Eq. 
(1.1) in its obvious time-dependent extension. This re- 
sponse, when looked at from the point of view of the 
system in its ground state, implies a redistribution of 
charge (wave functions), in other words, depletion of 
charge (holes) or accumulation of charge (electrons) in 
some places. This creation of electron-hole pairs can be 
seen explicitly in the corresponding equations. In fact, in 
direct (r) space, the dielectric function for an electron 
system can be written as 

e(r,r',w) = <?(r-r')- j" dr"v(r-r")P(r" y ,co) (1.3) 

where v is the bare Coulomb interaction and where a 
polarizability operator P(r",r',w) has been introduced. 
When P is zero, the system is not polarizable and hence 
the total potential is equal to the external one. Other- 
wise, P is, in general and on average, negative, i.e., act- 
ing against the external potential. The simplest approxi- 
mation for P is the independent (quasi)particle form: 



F iei .(r,r', W )==2 ifi-fj) 



<AKr)^*(r)(A,(r')^r(r') 
w— u>[j + irj 



To be precise, this calculation also includes a vertex correc- 
tion beyond the usual GW. 



(1.4) 

Here «,-,■= (e,— e ; ), /,■ are Fermi occupation numbers, 
and (/,/) label the states of energy e,- and e, obtained 
from some (for the moment, unspecified) equation for 
one-particle states. The small imaginary number irj 
leads to an imaginary part of P which is proportional to 
8(ej- €j- to); in other words, one can see the energy 
conservation for a photon to promoting an electron from 
state i to state /. 

The approximation for the polarizability P\qp in Eq. 
(1.4) as a sum over independent transitions has the form 
of the random-phase approximation (RPA; Adler, 1962; 
Wiser, 1963). In fact, the RPA was originally meant to 
describe a calculation performed within the (linearized) 
time-dependent Hartree approach, or Lindhard approxi- 
mation (Lindhard 1954), for a homogeneous electron 
gas. Later, the time-dependent Hartree approach was 
shown by Ehrenreich and Cohen (1959) to be equivalent 
to the diagrammatic bubble expansion for the dielectric 
function in many-body perturbation theory (see also 
Pines, 1963, and Pines and Nozieres, 1989). Here, the 
term RPA form will in the following indicate that the 
approximation (1.4) has been used for P. The "addi- 
tional" holes i and electrons can of course be calcu- 
lated within the framework of the one-particle excita- 
tions outlined above. It turns out that the dielectric 
function evaluated using Kohn-Sham orbitals if/j and ei- 
genvalues e t from LDA in Eq. (1.4) yields absorption 
spectra that are in quantitative, and sometimes even in 
qualitative, disagreement with experiments (Cohen and 
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Chelikowsky, 1988). Since the factor 5(e,— e t — w) im- 
plies that absorption occurs at eigenvalue differences, 
this has been in part attributed to the fact that the Kohn- 
Sham eigenvalues underestimate the quasiparticle gap, 
which leads to a redshift of the theoretical spectrum with 
respect to the experimental one. In fact, for example, in 
a small sodium cluster Na 4 the first allowed transition 
between Kohn-Sham orbitals occurs at 1.1 eV (Onida 
etal, 1995), whereas the first experimental absorption 
peak is seen at 1.8 eV (Wang et al, 1990). 

However, even when the quasiparticle gap problem is 
corrected by replacing the Kohn-Sham eigenvalues by 
quasiparticle energies calculated, for example, in the 
GW approximation (this approach will be named 
GW -RPA in the following), important discrepancies 
with experiment are found. In general, the redshift prob- 
lem becomes a blueshift problem — in the case of Na 4 , 
the first main absorption peak is shifted to 3.3 eV. This 
problem is not due to a failure of the GW approach. In 
fact, for Na 4 the ionization potential, i.e., the creation of 
a hole, is correctly described within GW (Reining et al, 
2000). In the latter case, however, one considers the ad- 
dition of a hole with respect to the ground state and not, 
as in the case of absorption, the addition of a hole and 
then the addition of an electron in the presence of that 
hole. An additional correction term to the potential is 
hence needed in the case of absorption, expressed by the 
change of the potential upon presence of the hole, and a 
self-consistent treatment of hole and electron altogether. 
One has therefore to expect an effective two-particle 
equation for the response function, including electron- 
hole interaction effects. This equation — the Bethe- 
Salpeter equation — is actually found via the rigorous 
treatment of the problem within Green 's-function theory 
(see Sec. IV.B). It requires the correct calculation of the 
quasielectron and the quasihole (for example, via GW) 
and contains an interaction term that mixes the formerly 
independent transitions. The electron-hole interaction 
hence modifies the expression for the polarizability in 
Eq. (1.4). These excitonic effects have been well known 
for some time, starting from a description based on a 
superposition of atomic excitations (Frenkel, 1931a, 
1931b; Peierls, 1932) and the insertion of the concept of 
"excitons" into the band picture (Wannier, 1937), and 
extending to detailed work in the 1950s like that of 
Heller and Marcus (1951), or of Haken (1956, 1957). An 
overview of the intense activity in the field at that time 
can be found in the book of Knox (1963). 

Sham and Rice (1966) established the link between 
exciton models and many-body theory, by deriving the 
effective mass approximation from the Bethe-Salpeter 
equation. The first realistic nonmodel calculation for va- 
lence excitons in a solid based on the Bethe-Salpeter 
equation was performed by Hanke and Sham, first using 
the time-dependent Hartree-Fock approximation (which 
is equivalent to introducing an unscreened electron-hole 
interaction; Hanke and Sham 1974, 1975) and later also 
introducing a static electron-hole screening (Hanke and 
Sham, 1980). Their calculation, based on the linear com- 
bination of atomic orbitals, managed to explain the 



qualitative features of the line shape of the spectra of 
diamond and silicon. The Na 4 cluster was the first real- 
istic system to be treated in the ab initio scheme starting 
from DFT and going through GW (Onida etal, 1995), 
with a particularly strong effect of the electron-hole in- 
teraction due to the finite size of the cluster: the first 
absorption peak was moved back by 1.5 eV from its GW 
position, with a result close to the experimental one. An 
increasing number of such ab initio calculations of exci- 
tonic effects in finite and extended systems (see, for ex- 
ample, Albrecht etal, 1997, 1998a, 1998b; Benedict 
etal, 1998a, 1998b; Rohlfing and Louie, 1998a, 1998b) 
have appeared since then, generally showing a big im- 
provement over the results of RPA and GW -RPA cal- 
culations. 

Of course these two-particle electron-hole calcula- 
tions are relatively cumbersome, and it has always been 
tempting to stay within the framework of DFT. One 
could in fact consider the DFT exchange-correlation po- 
tential V xc as an approximation to the self-energy and, 
knowing that self-energy and electron-hole interaction 
effects partially cancel each other, 4 hope that V xc to- 
gether with its variation with density (the functional de- 
rivative SV xc /Sp, i.e., the so-called exchange-correlation 
kernel f xc ), could yield improved optical spectra. This 
resembles the theoretically more rigorous TDDFT 
(Runge and Gross, 1984; Gross and Kohn, 1985), where 
the response of the electrons to a time-dependent exter- 
nal potential is derived by searching the extrema of the 
quantum-mechanical action functional, which leads to a 
time-dependent Kohn-Sham equation. It has turned out 
that TDDFT in the adiabatic local-density approxima- 
tion (TDLDA) often yields good spectra for finite sys- 
tems. For Na 4 and other sodium clusters, Vasiliev et al. 
(1999) have calculated absorption spectra in excellent 
agreement with experiment, i.e., reproducing peak posi- 
tions and relative heights within —10% (0.2 eV). This is 
not in contradiction with the fact that, as mentioned 
above, the Kohn-Sham eigenvalue differences are 
smaller than the absorption energies: TDDFT directly 
describes the evolution of the density under the influ- 
ence of a time-dependent external potential, which 
means that both the potential (yielding the Kohn-Sham 
eigenvalues) and variations of the potential (which can 
shift the excitation energies) are taken into account. 
TDLDA is also successful in describing electron-energy- 
loss spectra for bulk metals and semiconductors, but im- 
proves only very slightly the absorption spectra of solids 
with respect to RPA calculations (Gavrilenko and Bech- 
stedt, 1996). Therefore, in the field of TDDFT the main 
effort is directed toward finding better approximations 
for both the exchange-correlation potential and its varia- 
tion with the density (i.e., f xc ). To this end, one has to 
understand the main differences between finite and infi- 
nite systems, as well as between different spectroscopies. 



4 See the example of Na 4 in Onida et al. (1995), and the re- 
sults of the jellium model for metal clusters by Pacheco and 
Ekardt (1997; also Ekardt and Pacheco, 1995). 
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In fact, for example, absorption and electron- energy- 
loss spectra are derived in different ways from the same 
dielectric function. 

In the following we shall therefore concentrate in 
more detail on the dielectric function and its link to 
spectroscopy. In particular, we shall outline those fea- 
tures that are common to both TDDFT and Bethe- 
Salpeter approaches, before treating the two approaches 
separately and finally comparing them. 



II. TDDFT AND GREEN'S-FUNCTION APPROACHES: 
COMMON INGREDIENTS 

TDDFT is an extension of static ground-state density- 
functional theory. Also, Green's-function calculations of- 
ten start from DFT results. In fact, concerning the 
ground state, calculations based on density-functional 
theory using simple density functionals predominate to- 
day. Current DFT results have often surpassed those 
from standard ab initio quantum chemistry techniques 
like Hartree-Fock, configuration interaction, etc., which 
may require heavy compromises in their technical real- 
ization for systems that are not very small (e.g., full 
configuration-interaction calculations can hardly be 
done for more than five electrons; Jensen, 1999). De- 
spite the failures of DFT due to the approximations 
made for the exchange-correlation contributions, its use 
continues to increase because of its wide applicability 
and its favorable scaling with the number of atoms. 

In the following, we briefly review the basic ingredi- 
ents of DFT in view of its use as a starting point for 
spectroscopy calculations. For a more detailed descrip- 
tion we refer the reader to any of the numerous compi- 
lations published in recent years, e.g., Lundqvist and 
March (1983); Dreizler and Gross (1990); Seminario 
(1996). 

A. Density-functional theory for the ground state 

The ground-state energy of a system of interacting 
electrons in an external potential can be written as a 
functional of the ground-state electronic density. Com- 
pared to conventional quantum-chemistry methods this 
approach is particularly appealing, since it does not rely 
on a complete knowledge of the TV-electron wave func- 
tion but only on the electronic density. Of course, al- 
though the theory is exact, the energy functional is un- 
known and has to be approximated in practical 
implementations. Today's DFT starts with the theorems 
of Hohenberg and Kohn (1964) for a search of the elec- 
tronic ground state of an isolated system of N interact- 
ing electrons in an external potential V ext {r). The first 
theorem ("the density as the basic variable in the elec- 
tronic problem") establishes that the external potential 
V ext (r) is a functional of the charge density p(r), within 
an additive constant. The second theorem establishes 
the "energy variational principle for the density." These 
two theorems show that the problem of solving the 
many-body Schrodinger equation for the ground state 



can be exactly recast into the variational problem of 
minimizing the Hohenberg-Kohn functional with respect 
to the charge density. 

In practice, the available approximations for kinetic 
energy and exchange-correlation density-based func- 
tionals give moderate quantitative agreement with ex- 
perimental data (Dreizler and Gross, 1990). A much im- 
proved strategy has been presented by Kohn and Sham 
(1965) using orbital variables. The fundamental assump- 
tion is to introduce a reference system of noninteracting 
electrons in an external potential Vf f f (r) such that the 
ground-state charge density of the physical system coin- 
cides with that of the reference system. 5 With the intro- 
duction of the noninteracting system, the variational 
problem on p(r) is thus finally reformulated in terms of 
the following set of self-consistent Kohn-Sham equa- 
tions: 



1 



V 2 +V H (r) + V xc (r) + V (r) 



iA,(r) = e,iA,(r) 



(2.1) 



2 . Thus 



with the density given by p(r) = 2 1 = 1 |(/f,(r 
= V H + V xc +V , where V H is the Hartree potential, 
V X c( r ) = 8E xc [p]l Sp(r) is the exchange-correlation po- 
tential that contains all many-body effects, and V is an 
external potential stemming, for example, from the 
electron-ion interaction. 

The simplest approximation (still widely used) to E xc 
is the local-density approximation (Kohn and Sham, 
1965). The approximation is based on using the 
exchange-correlation energy density e'!" m of the homo- 
geneous electron gas (Ceperley and Alder, 1980; Ortiz 



etal, 1999), namely, u/i [p] = J e" x ° c m (p(r))p(r)dr. 
Hence one replaces the inhomogeneous electron system 
at each point r by a homogeneous electron gas having 
the density of the inhomogeneous system at r. The ra- 
tionale for this approximation is in the limit of slowly 
varying density. Unexpectedly, the domain of applicabil- 
ity of the LDA has been found to go much beyond the 
nearly free-electron gas and accurate results can be ob- 
tained for very inhomogeneous systems. Improvements 
over the LDA have been found, e.g., by the generalized- 
gradient approximations, in which the exchange- 
correlation energy density is a function not only of the 
electron density, but also of its gradient (see, e.g., Per- 
dew etal, 1996). Further improvements over the 
generalized-gradient functional should go in the direc- 
tion of nonlocal functionals, in order to describe better 
the inhomogeneity of the exchange-correlation hole, as 
in the case of the exact exchange potentials (Gross et ah, 
1996). 

The issue of calculating electronic excitations, how- 
ever, goes beyond the problem of finding a good ap- 
proximation to the ground-state exchange-correlation 



This definition is meaningful if the charge density is nonin- 
teracting v representable, so that one has a unique definition 
of the total energy Kohn-Sham functional. 
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functional. The problem is intimately linked to the fact 
that the one-particle eigenvalues in the Kohn-Sham 
theory have often been used, as well, to discuss the ex- 
citation spectra of solids, molecules and atoms. Al- 
though there is no rigorous justification for a direct iden- 
tification of Kohn-Sham eigenvalues with quasiparticle 
energies, the formal resemblance between Kohn-Sham 
and quasiparticle equations, such as the COHSEX 
equation [Eq. (1.2)], where 2(r,r',w) replaces <5(r 
— r')V xc (r) in Eq. (2.1) explains why the Kohn-Sham 
equation can be considered at least as a starting point. 

As pointed out above, once some one-particle-like 
theory such as Kohn-Sham has been established, re- 
sponse functions and spectra can be constructed, for ex- 
ample in the approximation (1.4). This is the subject of 
the next section. 

B. The inverse dielectric function and local-field effects 

As discussed above, information about the electrons 
and holes, and about their interaction, is contained in 
the dielectric function via the polarizability P of Eq. 
(1.3). However, further considerations are necessary in 
order to obtain spectra. In fact, independently of how P 
and e are calculated, the latter usually has to be in- 
verted, since one is concerned with the total potential 
Vtot =e ~ l Vext f° r a given external potential V ext , and 
not vice versa. Equivalently, the induced charge is given 
by Pind = XYext > where the reducible polarizability x an d 
e ~ 1 are linked by 

8~\t,t',ta)=S(x-x') + J dr"v(r-r") x (r" ,r' ,w) . (2.2) 

In the following, it is useful to describe explicitly two 
types of spectra containing neutral excitations, namely, 
photon absorption and electron-energy-loss spectra 
(EELS). In electron-energy-loss experiments, an elec- 
tron impinges on the sample and loses energy by excit- 
ing electron-hole pairs, plasmons, and other high-order 
multipair excitations. This energy loss is given by the 
imaginary part of the integral of the potential created by 
the electron, and the induced charge. When the poten- 
tial V ext due to an electron is taken proportional to a 
plane wave, this leads, for a momentum transfer q, to 
the loss function 

L(w)°c-Im| J dtdr' e~ iv x{**' ,<>>)e iv '\, (2-3) 

which is hence determined by the Fourier transform of 
the inverse dielectric function — Im[8 _1 (q,q,co)]. 

In a finite system, x a l so yields the photoabsorption 
cross section a, via 

er(w)= Imo(w), (2.4) 

where c is the velocity of light and Im a{a>) is the imagi- 
nary part of the dynamical polarizability, 

«(«)=- drdi'V ext (l,(o)x(l,l',(o)V ext (i',co). (2.5) 



In particular, for a dipolar external field V ext along the z 
direction, the corresponding component of the a tensor 
reads 

4 770) f 

(T zz {a>) = — Im I drdr' zx(x,r' ,w)z' , (2.6) 

which can also be obtained from the long-wavelength 
limit of e -1 : o-(<w) = lim q ^ (&>/c) Im[u(q)^(q,q,w)] 
along the field direction. The other components of the 
tensor can be similarly obtained. Note that in molecular- 
beam experiments, where the molecules are randomly 
oriented, one measures the trace of the tensor, that is, 

(r e x P (<>>) = l2i=(x, y ,z) or ii( (a )- 
In the solid, one has to distinguish between the loss 

spectra and the absorption spectra, which are instead 
described by the macroscopic dielectric function e M . Its 
relation to the microscopic dielectric function e of peri- 
odic crystals has been discussed by Adler (1962) and 
Wiser (1963) as well as by Ehrenreich (1966): 

1 

e M (&)) = hm— -, (2.7) 

q^O s G=0,G' = ( - < l' a} -' 

where 8 GG -(q,w):=8(q+G,q + G',w) is the Fourier 
transform to reciprocal space of e(r,r'), G is a recipro- 
cal lattice vector, and q belongs to the first Brillouin 
zone. The optical absorption spectrum is then given by 
the imaginary part e 2 (w) of e M (w). The dielectric con- 
stant 8 is the value of e M (co) at w = 0. If e is diagonal in 
G,G', s M is just e M = lim q ^ e G=0G , = , i.e., the spatial 
average of the microscopic dielectric function. This is in 
fact the case in the homogeneous electron gas. Other- 
wise, the microscopic dielectric function e(r,r') will de- 
pend explicitly on the positions r and r', and not simply 
on the distance |r— r'|. This is translated into the fact 
that the dielectric matrix in reciprocal space is not diag- 
onal. The fact that all elements of the matrix contribute 
to one element of its inverse reflects the so-called local 
field effects: these effects arise whenever the system un- 
der study is nonhomogeneous on the microscopic scale. 
In this case, for example, an external spatially constant 
perturbing field will induce fluctuations on the scale of 
the interatomic distances in the material, giving rise to 
additional internal microscopic fields. These concepts 
will be helpful in the next section [see the discussion 
after Eq. (2.23)], and in Sec. VI, where TDDFT and 
Bethe-Salpeter methods are compared, outlining the 
equations and the part of the kernel which are common 
to both approaches (see in particular Sec. VI.D.l). It is 
clear that the local field effects should in principle be 
included in both absorption and loss spectra, since the 
crucial quantity which appears is always the inverse di- 
electric function. 

Local field effects can shift peak positions, and they 
can be very important, as will be discussed later. How- 
ever, even when local field effects are taken into ac- 
count, the above conclusions about the need to go be- 
yond RPA or GW -RPA in the calculation of P, and 
hence e, remain valid (for example, the inclusion of local 
field effects does not generally lead to a significant im- 
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provement of the absorption spectra of solids; Louie 
et at, 1975). It should in fact be stressed that the concept 
of local field effects is independent of the approach used 
to calculate P and e. It is therefore useful to summarize 
here the notations that will be used throughout this re- 
view in order to characterize the different levels of ap- 
proximation: First, RPA and GW -RPA refer to a calcu- 
lation in which P has been approximated by the form 
(1.4), using DFT and GW eigenvalues, respectively, but 
the spectrum has then correctly been constructed from 
e _1 . When local field effects are neglected in the inver- 
sion, in the dipole limit Eq. (1.4) leads to the well-known 
Ehrenreich and Cohen (1959) expression for the imagi- 
nary part of the macroscopic dielectric function in terms 
of the matrix elements of the velocity operator v be- 
tween valence and conduction states: 6 

Im{e M ( w )}=- r X <<fc,M^) 2 5(e c -e B -a>). (2-8) 

W v,c 

Independently of the quality of the states if/ entering Eq. 
(1.4), the electronic excitations described by Eq. (2.8) 
are restricted to the generation of noninteracting 
particle-hole pairs. In the following this approximation 
will be called an independent-particle-random-phase ap- 
proximation macroscopic dielectric function (Del Sole 
and Girlanda, 1993). 

It is also interesting to point out that, on the other 
hand, if the matrix inversion is properly taken into ac- 
count, the formerly independent transitions mix. In 
other words, even if no electron-hole interaction is in- 
cluded in P, there is an effective electron-hole interac- 
tion showing up in s M . This is actually an electron-hole 
exchange term, as will be discussed later. 

The inclusion of local field effects is, in principle, 
straightforward: Inverting e and placing the Fourier 
transform of the inverse dielectric function in Eq. (2.7) 
directly yields the macroscopic dielectric function for 
each frequency co. For our purpose, however, it is con- 
venient to use a different formulation for the macro- 
scopic dielectric function, which is suitable for a subse- 
quent inclusion of excitonic effects. Moreover, and more 
importantly here, it will enable us to compare in a con- 
cise way (i) electron-energy-loss and absorption spectra 
and, later, (ii) TDDFT and the Bethe-Salpeter approach. 
In fact, one can show (see Hanke, 1978 and Appendix 
B.l) that e M can be constructed from a modified re- 
sponse function P, 

e M ( oo) = 1 - lim[>(q)oi 5 G=G' = o(q,<»)], ( 2 -9) 

where the matrix P GG< satisfies the Dyson-like screen- 
ing equation 



P = P + PvP. 



(2.10) 



6 This is none other than Fermi's golden rule in the dipole 
approximation and the one-electron picture for absorption 
processes (Pines, 1963; Fetter and Walecka, 1971). 



Here we have defined the modified Coulomb interaction 
u(q) G as follows: 

if G = 0) 

4tt 



u(q)< 



\q+G\ 



else 



In other words, the difference between v and v is "only" 
in the G=0 component. 7 Because of the long range of 
the Coulomb interaction, this term is divergent for van- 
ishing q, which explains its importance, as will be dis- 
cussed below, v, on the other hand, does not diverge and 
plays the role of a correction term, namely, the local 
field effects. In fact, from Eq. (2.10) it is clear that put- 
ting v to zero is equivalent to neglecting those effects. 

It is interesting to compare Eq. (2.10) to the corre- 
sponding Dyson-like equation for ^, which follows from 
Eqs. (2.2) and (1.3): 



x = P + Pv X , 



(2.11) 



the only difference is in fact the G = part of the bare 
Coulomb interaction. This remark turns out to be ex- 
tremely important, as will be discussed later in the com- 
parison of TDDFT and Bethe-Salpeter results for finite 
and infinite systems. One can put the common term v 
into evidence by writing Eq. (2.11) as 



XG,G'~ "g,G' 

+ 2 Pg,G"VG"XG",G' + ^ > G,O l; oA'0,G' 
G" 



(2.12) 



One can then immediately write down Eq. (2.12) and 
(2.10) for the case when local field effects are neglected, 
and take the macroscopic limit; two quantities Poo=Poo 
and A / oo = ^'oo/(l — ^"oo^o) are obtained, which describe, 
respectively, the absorption and the electron-energy-loss 
spectra of the system. Also in this simplified case it is 
clear that \ an d P are fundamentally different (x is 
screened, but P is not), and that this fact is entirely due 
to the seemingly tiny difference of the kernel of the 
Dyson-like equations, i.e., the difference between v and 
v (see detailed discussions and applications in Sec. 
VI.D). 

Until now, all considerations have been general, i.e., P 
has not been specified. Only when talking explicitly 
about P will there be a difference between TDDFT and 
the Bethe-Salpeter approach. However, in spite of — or 
because of — the differences that will arise, it is worth- 
while to anticipate a qualitative discussion of the struc- 
ture of the equations leading to the determination of P. 
This allows one to highlight the importance and conse- 



7 In a finite system, one can still have a wave vector k as a sum 
of a reciprocal lattice vector G and a vector q lying inside the 
first Brillouin zone, by creating a periodic array of "unit cells" 
of vacuum, each of them containing a copy of the system. The 
limit of infinite volume (isolated systems) for those supercells 
corresponds to the zero-volume limit for the first Brillouin 
zone, i.e., to an infinitely dense G-space. 
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quences of induced local and nonlocal potentials, with- 
out masking simple concepts by complicated many-body 
formulations. For the case of Hartree-Fock, many de- 
tails can be found in the reviews of McLachlan and Ball 
(1964) and Langhoff et al. (1972). 

C. Four-point kernels and effective two-particle equations 

An effective self-consistent Hamiltonian, which is 
made time dependent through the time dependence of 
the wave functions, in the presence of a time-dependent 
external potential, allows one to calculate the induced 
density matrix to first order (Hedin and Lundqvist, 
1969): 

p, w (r,r',w)= 2 ^„(r) 

n,n ' 

Jf n -f n ,)(n\H^(co)\n') , 



(2.13) 

where e n and if/ n {r) = {t\n) are the eigenvalues and 
eigenfunctions of the unperturbed Hamiltonian. Then x 
is determined through 



Pllul (r,r)= j dr' X (r,r')V ext (r'). 



(2.14) 



The perturbation H^ 1 ' not only is given by the external 
potential V ext , but also contains self-consistently the in- 
duced potential, which is proportional to the induced 
density matrix (schematically V ind = p ind SV ea l Sp). 

A typical effective Hamiltonian will have an effective 
potential V eii (r l ,r 2 ,t) that has some contributions that 
are local and depend on the density, i.e., 

V'/ oc {[p(iT,F)],r 1 ,r}<S(r 1 — r 2 ) (like the Hartree potential 
or the DFT exchange-correlation potential), and others 
that are nonlocal but directly proportional to the density 
matrix (like the Hartree-Fock exchange potential): 
p(r 1 ,r2,f) v v(r 1 ,r 2 ), where w(r 1 ,r 2 ) is some generalized 
interaction. For simplicity here we do not consider more 
complicated time-dependent interactions, but include 
the possibility of memory effects in V tuc , i.e., V[ oc can 
be a functional of the density at all (past) times 1. In this 
case, the induced potential is 



5(r 1 -r 2 )<5(r 3 -r 4 ) 



Vind( T i>*2>t)- dr 3 dr 4 dt 



X- 



SV loc ([p(riJ)],r u t) 
Sp(r 3 ,r 3 ,t') 



+ S(r l -r 3 )S(r 2 -r 4 )w(r l ,r 2 )S(t-t') 

X P(W( r 3> r 4> f ') 

dr 3 dr 4 dt' K(ri ,r 2 ,r 3 ,r 4 ,r,f') 

X Pm(/( r 3. r 4> f ')> 



which defines the four-point kernel K. Equations (2.13) 
and (2.15) can be put into a closed form for the induced 
density matrix, by introducing the four-point function 
4 Xo [see Eq. (B16)] such that Eq. (2.13) becomes 



Pmrf( r i> r 2,«) = J dr 3 dr 4 dr 5 dr 6 

x(l- 4 ^)- 1 (r 1 ,r 2 ,r 3 ,r 4 , w ) 

xV ext (r 5 ,w)S(r 5 -r 6 ). (2.16) 
Hence from Eq. (2.14) 

X(*i,*2,<*) = ) dr 3 dr 4 {l- 4 x Ky 1 (r 1 ,r 1 ,r 3 ,r 4 ,co) 



x Vo(r3,r 4 ,r 2 ,r 2 ,a>). 



(2.17) 



An equation of this kind is Eq. (4.9), derived below for 
the particle-hole response function in the Bethe- 
Salpeter scheme. Another one is the TDDFT equation. 
Note that in this case, since there are only local poten- 
tials, the S functions in the kernel (2.15) allow us to 
contract all equations immediately and to work with 
two-point functions from the beginning (see, e.g., 
Bertsch and Tsai, 1975; Bertsch and Broglia, 1994). In 
particular, the four-point function 4 xo can immediately 
be replaced by the two-point function P/gp(r 1 ,r 2 ) 
= 4 ^ (r 1 ,r x ,r 2 ,r 2 ) defined in Eq. (1.4), as in the main 
equation of TDDFT [Eq. (5.11)] in Sec. V.B. A compari- 
son of Eq. (2.17) and the four-point generalization of 
Eq. (2.11) allows one to get the four-point polarizability 
4 P from the equation 

(2.18) 



with 



4 /(*l ,*2 > r 3 M) = K ( r i ,*2 > r 3 M) - 4 v(ri ,r 2 ,r 3 ,r 4 ), 

and 8 

A v (r x ,r 2 ,r 3 ,r 4 ) = <5(rj - r 2 ) <S(r 3 - r A )v(r l ,r 3 ) . 

The two-point P is then 

P(r 1 ,r 2 ) = 4 J P(r 1 ,r 1 ,r 2 ,r 2 ). (2.19) 

Using the Green's-function formalism which will be in- 
troduced below (Sec. Ill), xo an d 4 Xo wm he expressed, 
respectively, as — i'G (12)G (21) and 

-iG (13)G (42). 

In general, V eff will of course contain the Hartree po- 
tential, and hence K will contain the variation of the 
Hartree potential with respect to the density, 8V H I Sp, 
which is just the contribution v appearing in the Dyson- 
like equation (2.11). This allows one to see the relation 
between the local field effects (see Sec. II. B) and the 
density variation of the Hartree potential: When local 
field effects are neglected, all spatial frequencies of the 



The definition of 4 y, which is used in the four-point equa- 



(2.15) tions involving P, is analogous to that of 4 v. 
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induced Hartree potential are neglected in the Dyson- 
like equation apart from the spatial frequency given by 
the external potential, i.e., apart from the G=0 contri- 
bution when the external potential is macroscopic. 9 In 
finite systems, macroscopic potentials do not play a role, 
and the distinction is meaningless — neglecting local field 
effects becomes equivalent to completely neglecting the 
induced Hartree potential. Therefore Eq. (2.6) is at the 
same time proportional to the macroscopic loss function 
and the to limit of Eq. (2.9) for the case of an infinitely 
dense G space. 

To summarize, one will in general find a four-point 
equation for S = 4 x or S = 4 P: 



S = 4 Xo + 4 XoKS, 



(2.20) 



where K contains (K = v+f) or not (K—v+f) the 
long-range term u(G=0), depending on whether S = 4 x 
or S = 4 P. Any density-matrix-dependent effective po- 
tential will lead to a four-point equation for S. In order 
to get the two-point functions, S has to be contracted 
like in the expression for the macroscopic dielectric 
function, given by 



e M (w) = l-lim 



q^O 



X 



v(q) j drdr'e-^ ^-"'^ 4 P(r,r,r' ,r' ; w) 



(2.21) 

If no nonlocal terms are present in V eff (e.g., when V e£f 
= V Hartree + V? C FT )' Eq. (2.20) can be contracted immedi- 
ately, and one can work with two-point functions only. 

Generally, in order to solve Eq. (2.20), one has to in- 
vert a four-point function for each frequency. This has 
been done even for solids, for example in the case of a 
Fock term in the effective Hamiltonian, in the early 
work of Hanke and Sham (1974, 1975). However, for a 
comparison between TDDFT and the Bethe-Salpeter 
equation method it is more convenient to present the 
equations in a different way, also often used. This 
scheme, described in the following, has the advantage of 
putting the two-particle nature of the problem into evi- 
dence. In fact, it has been shown (see Appendix B.2) 
that the four-point equation for S can be transformed to 
an eigenvalue problem involving the effective two- 
particle Hamiltonian: 

■:(e n -e ni )S nin3 S n2 „ 4 

+ (/n 1 _ /« 2 )-^(n 1 n 2 ),(« 3 « 4 ) • (2.22) 

The indices « ; refer to the fact that matrix elements 
have been taken with respect to four eigenfunctions of 
the starting effective static one-particle Hamiltonian. 
Hence H 2p is diagonalized, and from its eigenvalues E x 



M 2p 

(n 1 n 2 ),(n 3 n 4 ) 



9 The local field effects can also be included through a direct 
calculation of the variation of the Hartree potential in re- 
sponse to an electric field as done in density functional pertur- 
bation theory (Baroni et al., 2001 and references therein). 



and eigenstates A n ^ ni the four- and two-point quantities 
of interest are constructed. For example, the macro- 
scopic dielectric function is given by 



8 M (w) = l-limu(q)2 2 ,q 'V2> 

q^0 \«1.«2 



A 



(«1«2) 



w-E k + irj 



X 2j ( n 4\e 



iqr 



l«3Mw (/», 



-fn 4 ) 



(2.23) 



Here N k ,\' is an overlap matrix (see Appendix B.2). 
Since it can be shown that the n, always appear as pairs 
consisting of one occupied and one empty state (see Ap- 
pendix B.2), this reflects the very physical approach of 
rewriting the problem in terms of mixing between differ- 
ent independent-particle transitions between occupied 
and unoccupied eigenstates of a one -particle Hamil- 
tonian. Equation (2.23) should be compared to Eq. (2.8); 
the difference is caused by a change in excitation ener- 
gies [from (e c — e v ) to E k ], but also, very importantly, by 
the coefficients A k which mix the formerly independent 
transitions. As pointed out above, even if only the Har- 
tree potential is used in the derivation of K, i.e., if one 
considers only local field effects and no further 
exchange-correlation contributions, the electron-hole 
Hamiltonian is not diagonal, and the coefficients ^4 X will 
mix the transitions. This way of formulating the problem 
is hence also appealing because the knowledge of the A K 
allows one to interpret spectra in terms of mixing of 
transitions. 

One can thus define a general two-particle Hamil- 
tonian which (i) yields, depending on its interaction po- 
tential, the inverse dielectric function or the macro- 
scopic one and (ii) does so for systems that are described 
by an effective one-particle Hamiltonian either depend- 
ing on just the density, or containing more complicated 
functionals involving the density matrix. 

The corresponding equations have been derived start- 
ing from the time-dependent density matrix in Eq. 
(2.13). In order to talk now in detail about the possible 
effective one-particle Hamiltonians and derived kernels, 
one is naturally led to rigorously introduce in the follow- 
ing a sort of time-dependent density matrix for electrons 
and holes, which will allow the effective Hamiltonian 
and its derived response functions to be obtained on the 
same footing; one is led to introduce the Green's 
function. 10 



Note that the history of modern band-structure calculations 
started with the Hartree theory, where the essential ingredient 
is the electronic density, p(r) = 'Z s tfff (r)ip s (t). Exchange ef- 
fects (Hartree-Fock) are then included through the nonlocal 
density matrix p(t,r') = 'Z s if/f (r')i[r s (r), whereas to introduce a 
time dependence, the natural way is to do it in the phase of the 
wave functions, which leads to a time-dependent density ma- 
trix: p(r,r',f'-f) = 2 s i/f*(r')(/' J (r)e'' f * ( '~'' ) . This quantity is ac- 
tually a sort of one-hole Green's function, which will show up 
in the following in the more formal discussion of the problem. 
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III. GREEN'S-FUNCTION THEORY 

The Green's-function technique is indeed extremely 
useful: it allows one to tackle directly the problems of 
calculating excitation and ionization energies, ground- 
state energies, transition matrix elements, absorption co- 
efficients, and dynamical polarizabilities, as well as elas- 
tic and inelastic electron cross sections. Furthermore, 
self-consistent perturbation theories can be formulated 
in terms of the Green's function. The technique is fully 
discussed in many textbooks and we refer the reader to 
those for details (Abrikosov et al, 1963; Hedin and Lun- 
dqvist, 1969; Fetter and Walecka, 1971; Landau and Lif- 
schitz, 1980), but in the following we outline the equa- 
tions needed in the discussion. For further reading about 
the subtleties of the mathematical concepts involved, see 
the review by Farid (1999a) where this many-body ap- 
proach is developed in detail together with the approxi- 
mations to Hedin 's equations discussed below. 

A. The concept of Green's functions and the self-energy 

The one-electron Green's function G is defined (at 
zero temperature) as an expectation value with respect 
to the many-electron ground state \N), 

G(xt,x't') = -i(N\Ti/f(xt)^(x't')\N), (3.1) 

where <p(xt) is the field operator in the Heisenberg pic- 
ture, x stands for three space coordinates (r) plus one 
spin coordinate (£), and T is the time-ordering operator. 
In this equation, iffl(x,t)\N) represents an 
(N+l) -electron state in which an electron has been 
added to the system at point r and time t. When t'<t, 
the many-body Green's function gives the probability 
amplitude to detect an electron at point r and time t 
when an electron has been added to the system at point 
r' and time t' . When t'>t, the Green's function de- 
scribes the propagation of a many-body state in which 
one electron has been removed at point r and time t, 
that is, the propagation of a hole. 11 G is closely related 
to fundamental properties like charge and spin density 
or the total energy, and derived spectra like photoemis- 
sion or Compton scattering. Thus, for example, the 
charge density is obtained through 



p(rt)- 



i lim. 



G(xt,xt+r)d%. 



Similarly, the total electronic energy can be deduced 
through the Galitskii-Migdal formula (Galitskii and 
Migdal, 1958) 



E= — | dx lim,/^ f + 



dt 



- ih(x) 



G(x,t,x' ,t') x < 



(3.2) 



11 When time-reversal symmetry holds, the Green's function is 
symmetric with respect to the interchange of the spatial coor- 
dinates. This general statement is also true for the density ma- 
trix. 



where h(x) is the density-independent part of the effec- 
tive Hamiltonian. By inserting a complete set of (N 
+ 1)- or (N— 1) -particle states (depending on the time 
ordering) between the field operators in Eq. (3.1), one 
obtains the Lehmann representation of the Green's 
function in terms of the amplitudes, 

\N\ftx)\N+U), e s >fi 
(N-l,s\if/{x)\N), e s <fji 

and in terms of the real energies 

-E(N-l,s)]) for e s >pi (</*): 

fsix)ff(x') 



/,(*) = 



(3.3) 
■E(N+l,s) 



E(N) (or [E(N)- 
G(x,x' ,w) = 2 



a>-[e s + ir/sgn(/j,-e s )] ' 



(3.4) 



where the small imaginary part 97 is needed for the con- 
vergence of the Fourier transform to frequency space. 
The poles of G hence correspond to electron addition 
and removal energies. In particular, the so-called spec- 
tral function A(a>) :=|(l/7r)Im G(w)| becomes 



A(jc,x';w) = 2 f s (x)ff(x')S(co-s s ). 



(3.5) 



In the noninteracting case, where \N), \N+1) and |JV 
— 1) are simply Slater determinants, only states with 
|A f +l,s) = cJ|7V) (c^ is a creation operator) lead to non- 
vanishing Lehmann amplitudes. These amplitudes and 
the Lehmann energies are then equal to the eigenfunc- 
tions and eigenvalues of the corresponding one-electron 
Hamiltonian, and the spectral function consists of a set 
of S peaks at those eigenvalues. Hence each peak corre- 
sponds to a particle. Furthermore, it can be easily seen 
that the Green's function reduces to the time-dependent 
electron or hole density matrix introduced in Sec. II. C. 
When the electron-electron interaction is turned on it is 
no longer possible to write \N — 1) or |JV+1) simply in 
terms of pure one-electron Slater determinants. Instead, 
the many-body wave function is a linear combination of 
many of them. For example, the (N + 1) -particle wave 
function can be expressed as \N+l,s) = ^ lm a s lm c] n \Ni), 
where \N/) is an excited state of the iV-particle system. 
Hence there will now be more nonvanishing contribu- 
tions to the spectral function (3.5). If these contribu- 
tions, merged together, form a clearly identifiable main 
structure, which can be thought to derive from a S peak 
when the interactions are switched off, one can still 
work in a particle-like picture, associating each peak 
with a "quasiparticle." 12 With respect to an 
independent-particle 8 peak, for the quasiparticle one 



The concept of quasiparticles was introduced by Landau 
(Landau, 1957a, 1957b, 1959; Landau and Lifschitz, 1980) in 
the Fermi-liquid theory of ordinary metals as a one-to-one cor- 
respondence between low-energy excitations of a free Fermi 
gas and those of the interacting electron liquid. A quasiparticle 
can be considered as the combination of a real particle (elec- 
tron or hole) and a cloud of virtual electron-hole pairs sur- 
rounding it. Through mutual interactions quasiparticles can de- 
cay into other quasiparticles, leading to a finite lifetime. 
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finds a spreading of the oscillator strength. Therefore 
the broadening carries valuable information about 
many-body correlation effects in the interacting system. 

Another important quantity for the analysis of the op- 
tical spectra of interacting electron systems is the two- 
particle Green's function, defined as 

G(l,2,l',2')=-(N|r i A(l) ( A(2)(A t (2') ( A t (l')|iV), 

(3.6) 

where space, time, and spin coordinates are indicated in 
abbreviated form by numbers, i.e., 1 = (r 1 ,t\ ,£i). This 
function provides information about processes involving 
two-particle transitions and their interaction through the 
time ordering of the four times appearing in the equa- 
tion. For the purpose of this review, only a few specific 
time orderings need be considered, in particular for 

t\ ,t[>t 2 ,t 2 

G(l,2,l',2')=-2 X S (1,V)X S (2,2') 

s 

in terms of the hole-particle (Bethe-Salpeter) ampli- 
tudes X s (l,V) = (N\Ttfi(l)^(l')\N,s) and X S (1,V) 
= (N,s\Ti//(l)i//\l')\N). Also the case h,t[<t 2 ,t' 2 de- 
scribes particle-hole amplitudes, whereas other time or- 
derings will lead to particle-particle or hole-hole pair- 
ings. This point is relevant here because only for hole- 
particle pairings do the states belong to the iV-particle 
system. Therefore one can define a specific hole-particle 
Green's function G hp considering only those two time 
orderings (/ t ,t[^t 2 ,t' 2 ): 

G( 1,2,1 ' ,2' ) = G hp { 1,2,1 ' ,2' ) + other orderings. 

(3.7) 

The "other orderings" terms correspond to poles of the 
two-particle Green's function that differ from 
±[E(N,s) — E(N)] and that are important for different 
types of spectroscopies, e.g., Auger in the case of G hh 
(Cini, 1977, 1979). 



1 . Connection to spectroscopic measurements 

In the Introduction we sketched the main spectro- 
scopic measurements that can be described by the ap- 
proaches presented in this review. It is now useful to 
build a bridge between those qualitative considerations 
and the abstract mathematics of the previous subsection. 

In fact, the building block for the description of any 
spectrum involving the interaction of radiation and mat- 
ter is the probability V for an excitation from the initial 
state |^o) (typically the ground state) to a set of final 
states l^y). According to the measurement, the initial or 
final states contain an additional photon. This process is 
well described by Fermi's golden rule: 



terms in the vector potential A, A= (e/2mc) (Ap 
+ p.4). High photon fluences (multiphoton processes) 
are excluded here. 

From the transition probability, one can construct the 
spectra. In the case of photoemission (see Fig. 1), one 
measures the photocurrent, which is proportional to the 
probability per unit time of emitting a photoelectron of 
momentum k, when the target is irradiated with photons 
of frequency w. Hence one looks explicitly at the photo- 
electron. As a consequence, the simple approximations 
for photoemission decouple the photoelectron from the 
target, usually treating it as a free propagating electron 
state. One then finds that the photoemission spectrum 
carries directly accessible information about the target 
itself, namely, about its electronic structure, or to be pre- 
cise, the energy of holes (or electrons, for inverse pho- 
toemission). The final simple relation between the pho- 
toemission spectrum and the electronic structure of the 
sample is based on the classical three-step model (see, for 
example, Almbladh and Hedin, 1983, and Kevan, 1992): 
the photoelectron is excited by the photon, loses energy 
on its way to the surface, and finally has to pass the 
surface to propagate (as a free particle) to the detector. 
Only the first step, the excitation of the photoelectron, is 
contained in the so-called intrinsic spectrum, whereas 
losses on the way out of the target are called extrinsic. 13 

Starting from the transition probability [Eq. (3.8)], 
one can write down the photoelectron current of photo- 
emission, using the initial-state energy E =E N + to, i.e., 
the sum of the iV-particle ground-state energy and the 
photon energy. In the intrinsic approximation, the final 
state \k;N—l,f) is approximated by cl\N—l,f) with en- 
ergy €k+E N -i f, where k and e k are the momentum and 
the energy of the photoelectron. The photoelectron cur- 
rent reads then (Almbladh and Hedin, 1983) 



k 



X 



<W-l,/k k 2 Aijc} Cj \N) 



(3.9) 



Here e f =E N -E N 



and ~2 i jA i jcJcj is the photon op- 
erator in second quantization. Using the argument that 
the high-energy one-electron state k does not change the 
state of the (N — 1) -electron system, and that it makes a 
negligible contribution to the virtual one-electron exci- 
tations contained in the ground state N (that is Cy\N) 
= 0), one can write the photoelectron current as 



J k (w)- 



k 
4V 



2 A, 



;A ij (€ k -w)A 



jk j 



(3.10) 



7>=2t72 \(Vf\A\y )\ 2 S(E f -E ), 



(3.8) 



where A t j are the matrix elements of the spectral func- 
tion (3.5) between one-electron orbitals, 

Ajj(e k — w)= j dxdx't//f(x)if/j(x')A(x,x',e k —co). 

(3.11) 



where A describes the perturbation due to the photon 
field. For a one-photon process, neglecting quadratic 



13 For high-energy photoelectrons this contribution is, in gen- 
eral, much smaller than the intrinsic part. 
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If one further assumes that the diagonal elements of A 
are dominating (which is exact for independent elec- 
trons) and that the A are constant, the photoemission 
spectrum is proportional to the trace of the interacting 
spectral function (also neglecting electron-phonon cou- 
pling and impurity effects). 

From Eq. (3.5) one can see that the spectral function 
has a direct link to the density of available states for the 
addition and removal of electrons. In fact, for indepen- 
dent particles the trace of the spectral function is pro- 
portional to the density of occupied and empty states of 
the system, respectively, and provides a good approxi- 
mation for the direct and inverse photoemission spectra. 
As we have seen, the same is true for the interacting 
electron case: the trace of A(x,x'\w) can be accessed 
experimentally. In this way, the quasiparticle spectrum is 
directly measured in experiments such as direct, inverse 
and time-resolved two-photon photoemission (see Fig. 1 
for a schematic diagram of the first two). Still, final-state, 
extrinsic, and surface (transport) effects, neglected in 
this discussion, can play an important role in some spe- 
cific experimental configurations. 

The absorption coefficient, by contrast, is defined as 
the ratio of the energy removed from the incident beam 
per unit time and per unit volume to the incident flux. 
This energy is absorbed thanks to transitions from the 
initial to any final state, with the sole constraint of en- 
ergy conservation. There is, in general, no simple as- 
sumption concerning the final state as there is in photo- 
emission, and one has to deal with holes and electrons at 
the same time. It would not be a good approximation for 
the final state to decouple the excited electron from the 
others, nor to describe it by a scattering (plane-wave) 
state. One has therefore to sum over all possible 
jV-particle excited states in order to obtain the absorp- 
tion coefficient. The energy differences appearing in the 
<5 function are now E N t— (E N + co), and the sum over all 
states leads to a density-density correlation function (or 
density response function). 14 The latter can be obtained 
from the time-ordered (or retarded) two-particle 
Green's function G hp (l,2,V ,2') [Eq. (3.7)], by contract- 
ing some of its arguments and by taking suitable aver- 
ages or integrals of the resulting (microscopically vary- 
ing) functions. That is, the microscopic polarizability of 
Eq. (1.3) or the density response function entering in 
Eqs. (2.2) and (2.5) are obtained from a G hp where the 
limit (1')->(1 + ); (2')^(2 + ) has been taken. In the 



Indeed, using the standard Kubo formula for the density- 
density response function (Abrikosov et ai, 1963) ^(x,f;x',0) 
= -i8(t)(N\[p(x,t),p(x'fi)]\N), together with the clo- 
sure relation of TV-particle states \N,f) and the time evolu- 
tion of the density operator p(x,t) = e' Ho 'p(x,Q)e~' H °', 
one finds the response function ^(x,r;jc',0) 
= -id(t)[i: f B f (x)Bf(x')e i< ' E N- E N,f)'-c.c.l where B f (x) 
= {N\p(x,t)\N,f) and E N j is the energy of the /th N-particle 
excited state. Then its time Fourier transform has poles at the 
true excitation energies of the iV-particle system, a>=±(E N 



RPA approximation, P\qp is derived by taking G hp as 
the product of two one -particle Green's functions: 
G(1,2,1',2') = G(1,2')G(2,1'), which becomes 
G(1,2)G(2,1) after the contraction. More details are 
given in Sec. IV.B.2. 

2. Dyson's equation: the self-energy operator X 

The full one-particle Green's function cannot be com- 
puted exactly, and one therefore needs a method to de- 
rive a physically sound approximation. In the simple 
case of the Hartree-Fock approximation, G is written 
directly in terms of the Hartree-Fock eigenvalues and 
eigenfunctions as G HF (x,x' ,w) = ^L s ip s (x)ip*{x')l{w 
— [e s + irfSgTi(fi— e s )]}. The inverse of G HF is hence 
linked to the Hartree-Fock Hamiltonian H : 

( G HF ) ~ 1 (x ,x ' , w) = w - H HF (x ,x'). (3.12) 

For the general case, it is therefore reasonable to write 
the frequency Fourier transform of G - 1 as 

G~ 1 (x,x',(o) = {S(x — x')[co — H H (x)] — '%(x,x',(i>)}, 

(3.13) 

where H H contains the kinetic-energy operator, the ex- 
ternal and the Hartree potential, and 2 is supposed to 
include the Fock exchange and the rest of the Coulomb 
interaction (correlation). This means that in this Dyson 
equation for G the exchange (the Fock term) and cor- 
relation contributions are included in an operator that is 
in principle unknown, the self-energy X. 15 Since the 
right-hand side of Eq. (3.13) contains a>, i.e., the Fourier 
transform of a time derivative, implicit expressions for 2 
can be found through an evaluation of the equation of 
motion for the Green's function. This scheme gives rise 
to a set of coupled differential equations in which 
higher-order Green's functions appear. Those equations 
can be exactly decoupled in only a few cases, and one 
has to resort to approximations to get a practical, solv- 
able scheme. However one can sum up the series of 
coupled equations for the Green's functions in a pertur- 
bative way and recast all the information in the operator 
X, as we shall see later. 

Suppose for the moment that 2 is known. Expressing 
G in terms of the Lehmann amplitudes and energies 
[Eq. (3.3)], one finds, after taking (for a discrete level) 
the w^e s limit, 

J {S(x-x")[s s -H H (x)]-X(x,x",B s )}f s (x")dx"=0, 

(3.14) 

i.e., the Dyson equation for the quasiparticle energies e s 
and amplitudes f s . The amplitudes f s form a complete 
set but are nonorthogonal, since the self-energy opera- 
tor is energy dependent. This also implies that the equa- 
tion can have more solutions than in the static case. In- 



See Farid (2001) for a detailed discussion of the general 
properties of the dynamical self-energy operator, including an 
analysis of the shortcomings of the usually applied approaches. 
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stead, if the electron-electron interaction is smoothly 
switched off, those amplitudes tend to form the standard 
orthonormal set of an independent-particle system. This 
leads again to the situation described above, with the 
difference between the isolated independent-particle 
peaks and the broad quasiparticle ones containing a cer- 
tain number of nonvanishing Lehmann amplitudes. Dy- 
son's equation (3.14) indeed has a one-electron form 
(single-particle Schrodinger-like) but the operator X has 
very nontrivial correlation effects built in and is far from 
a mean-field approximation. 

In order to deal more easily with the pole structures 
of G, 2, etc., it is convenient to generalize all quantities 
from the real frequency axis w to the complex plane z, 
i.e., to make an analytic continuation. The one-particle 
Green's function is then written in terms of left and right 
solutions 



G(z) = X 



fl(z)^'U) 
z-E s (z) 



(3.15) 



defined by the general Dyson equation to be solved in 
the full complex plane, 

E s {z) = z, 
where 



[H +X(z)W s (z) = E s (z)cf> r s (z), 
[H + Mz)Vct>'*(z) = E?(z)ct> l *(z). 



(3.16) 



The self-energy 2 here plays the role of an effective 
potential for an electron or hole added to the system 
which is, in general, complex, nonlocal, and energy de- 
pendent. This potential arises from the exchange and 
from the response of the rest of the electrons to the 
presence of the additional particle. 2 ( u>) is real when w 
is set to an energy value below the first inelastic thresh- 
old of the system. 

Thus one has achieved an exact one-particle picture at 
the price of introducing a complicated effective poten- 
tial. The non-hermiticity and energy dependence of 2 
imply that the quasiparticle wave functions are energy 
dependent and not necessarily orthogonal; however, 
both the Green's function and Dyson's equations admit a 
biorthogonal spectral representation (in terms of left 
and right eigenfunctions). 16 The important point to note 
is that now one is dealing with the energies E s (z), which 
are in principle complex. 

The mathematical details of the analytic continuation 
are nontrivial and have important consequences. For a 
rigorous and complete treatment we refer the reader to 
the review paper of Farid (1999a). It is, however, useful 
to illustrate the connection between the Lehmann rep- 
resentation and the analytic continuation by an example: 



16 The biorthonormal and Lehmann representations are both 
exact. Even if they look similar they are not identical, since in 
the former case one works with the analytic continuations 
(Farid, 1999a). 



Suppose that some matrix element of a Green's func- 
tion is given in the Lehmann representation by 



G(«) = E 



e s + irj 



(3.17) 



For simplicity, let us consider a case in which g(s) 
= (1/tt)E 2 /[(s-E 1 ) 2 + eI] and e s =s. Given this par- 
ticular form of g(s), taking the thermodynamic limit 
(i.e., making s a continuous variable) yields a spectral 
function A{w) = 1Itt\ Im[G(w)]| with a Lorentzian form: 



A(w) = 



77 {(0-E X ) 2 +El' 

On the other hand, we can also write 



A(«) = 



1 



Im 



1 



a>-(E 1 + iE 2 ) 



(3.18) 



(3.19) 



thus the Lorentzian, which in the beginning was given by 
the imaginary part of a function containing a series of 
infinitely close-lying poles on the real axis (a branchcut), 
in this particular case can also be described by the imagi- 
nary part of another function, defined in the whole com- 
plex plane, having a single, but complex, pole (E 1 
+ iE 2 ). The position of the peak of A ( w) is given by the 
real part E\ , and the imaginary part E 2 gives its width. 
This simple example demonstrates the connection be- 
tween the Lehmann and the complex-pole representa- 
tion of the Green's function, in the case of a single qua- 
siparticle pole. It also explains why the complex 
representation is advantageous — in this example it re- 
places the search for a branchcut by the search for one 
complex pole. Moreover, it illustrates the meaning of E\ 
and E 2 : the real part of the quasiparticle energies solu- 
tion of Dyson's equation gives the band structure, and 
the imaginary part the quasiparticle damping (electron 
dynamics). The solution can be real valued only for ex- 
citations with infinite lifetime [Im2(z = e 5 )=0; the dis- 
crete part of the spectrum]. However, for energy ranges 
where e forms a continuum, 2(x,x';e) is also complex 
and Dyson's equation still has a solution with a complex 
z. In principle, this non-Hermitian eigenvalue problem 
can be solved for each and every value of z in the com- 
plex plane (with the corresponding set of eigenvalues 
and eigenvectors). Finally, this simple example is meant 
to explain why the existence of complex quasiparticle 
energies is not in contradiction with the fact that e s ap- 
pearing in the Lehmann representation of the Green's 
function are always real quantities. 

In order to extend the quasiparticle picture, one can 
return to the description of a photoelectron spectrum in 
terms of the Green's function and self-energy operator 
(in its analytic continuation, i.e., using complex ener- 
gies), by considering the matrix elements of Eq. (3.15): 



G „{(»)-- 



1 



(3.20) 



u> — Ei{a>) 

If one expands this expression for energies w close to the 
generally complex quasiparticle energy E t , one obtains 

G "(«)~ ~ n r ' 
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where Z, is the complex renormalization factor, 
1 



1 



do) 



(3.21) 



The «th matrix element of the spectral function for w 
close to Ei then reads 

1 llmf/ReZ, — (w — ReEjOImZ,) 

A u {w)= . 

7T (w-Re£ I ) 2 +(Im£ i ) 2 

If |Im£ ; | is small, A u {a>) has a sharp Fano peak at a> 
= ReEj of width T^Imi?,! and strength |ReZ,-|, with 
the asymmetry set by ImZ, . Z, is hence a measure of 
the strength of the quasiparticle pole and determines the 
specific shape of the corresponding excitation in the 
spectral function. Even for weakly correlated systems 
like sp metals and valence semiconductors | Re Z,| is far 
from 1 and can typically be in the range 0.6-0.9 (Hedin 
and Lundqvist, 1969; Hybertsen and Louie, 1986; Godby 
et al, 1988). The sum of all resonances, i.e., the trace of 
A, will finally determine the specific shape of the pho- 
toelectron spectrum. 



B. Hedin's equations 

As discussed in the Introduction (Sec. I.B), when go- 
ing beyond Hartree-Fock one has to deal with electron 
relaxation and correlation effects, which can also be de- 
scribed as dynamical screening effects. In fact, one can 
write an expansion of the self-energy in terms of a 
screened Coulomb potential W, by demanding that the 
total energy be stationary with respect to variations of 
the Green's function (Hedin, 1965). In this process the 
self-energy should be seen as a functional of the Green's 
function (2[G]). The use of W instead of the bare Cou- 
lomb potential is physically more sound: conventional 
many-body perturbation theory suffers in fact from vari- 
ous convergence problems, which have to be bypassed 
with care (see, e.g., Farid, 1999b). 

In fact, interactions in a real many-body system are 
screened to a large extent, and therefore W should be a 
much better behaved quantity for the development of a 
perturbation expansion than the bare Coulomb poten- 
tial. The latter is known to lead to convergence prob- 
lems in the range of densities of normal metals and 
semiconductors (Hedin, 1965, 1999; Hedin and Lun- 
dqvist, 1969, 1971). This expansion allows one to write 
the many-body problem as a closed set of five integral 
equations introduced by Hedin (1965), which relate the 
Green's function, self-energy, polarization propagator, 
and vertex function. The equations can be obtained by 
introducing a local time-dependent source term in the 
Hamiltonian which directly couples to the particle den- 
sity. This new source term is set to zero once the final 
equations are obtained. In this way Hedin derived a self- 
consistent scheme written in terms of the following set 
of coupled equations [here, as before, 1 = (r 1 ,t l and 
v stands for the bare Coulomb interaction]: 



2(12) = i J" G(13)r(324)W(41)d(34), (3.22) 
W(12) = u(12) + | w(13)P(34)W(42)d(34), (3.23) 



P(12) = -i G(13)G(41)r(342)d(34), 



(3.24) 



T( 123) =£(12) £(13) 



/ 



£2(12) 
£G(45) 



G(46)G(75)r(673)J(4567), 



(3.25) 

where P(12) is the time-ordered polarization operator, 
W(12) is the dynamical screened interaction, and 
T(123) is the vertex function. The fifth equation is the 
Dyson equation, which links G and £ [see Eq. (3.13)]. 
This way of writing the equations is in fact appealing, 
since it highlights the important physical ingredients: the 
polarization [Eq. (3.24)], which contains the response of 
the system to the additional particle or hole, is built up 
by the creation of pairs of particles and holes (the two 
Green's functions). The vertex function T contains the 
information that the hole and the electron interact. T, in 
turn, is determined by the change in the potential upon 
excitation [Eq. (3.25)]. 

C. The iterative approach 

Hedin's equations together with Dyson's equation 
form a set of equations that must in principle be solved 
self-consistently for G. This means that the Green's 
function used to calculate the self-energy should coin- 
cide with the Green's function obtained from the Dyson 
equation with the very same self-energy. It is obvious 
that this is a very difficult task, and that one must in 
practice find a simplification of Hedin's equations. Since 
it is not possible to find a straightforward, well-defined 
and convergent perturbation expansion in some small 
parameter, the approach to simplifying the equations is 
somewhat arbitrary and needs to be analyzed in detail. 
A self-consistent calculation of the interaction of low- 
energy electrons with an electron gas was first carried 
out by Quinn and Ferrell (1958). They performed a self- 
energy calculation of electron-electron scattering rates 
near the Fermi surface and derived a formula for the 
inelastic lifetime of hot electrons which is exact in the 
high-density limit. These free-electron-gas calculations 
were extended by Ritchie (1959) 17 and Quinn (1962) to 
include, within the first Born and random-phase ap- 
proximations, energies away from the Fermi surface, 
and by Adler (1963) and Quinn (1963) to take into ac- 
count the effects of the presence of a periodic lattice 
and, in particular, the effect of virtual interband transi- 



17 The 1/2 factor in front of z 2 in the expansion of f x just 
before Eq. (6.15) of this reference must be replaced by 1/3, as 
done in a subsequent paper (Ritchie and Ashley, 1965). 
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tions (Quinn, 1963). In their original work, Quinn and 
Ferrell (1958) had already given the GW expression for 
2, but did not recognize its usefulness beyond the case 
of an electron gas and only evaluated the results as a 
density (r s ) expansion. Further studies for quasiparticle 
excitations and lifetimes in a homogeneous electron gas 
were done by Hedin (1965) and Lundqvist (1967, 1968). 

1. First iteration step: the GW approximation 

The simplest approximation to T(123) assumes this 
operator to be diagonal in space and time coordinates, 
r ( 123) = <5(12) 8(13): the so-called vertex corrections 
are neglected. This implies that the irreducible polariz- 
ability is given by noninteracting quasielectron- 
quasihole pairs, 

P(12)=-/G(12)G(21 + ). (3.26) 

It can be evaluated explicitly by using Eq. (3.4). After a 
Fourier transformation, and taking into account the 
small imaginary parts of the energies, one obtains the 
formula (1.4) in its time-ordered version; this approxi- 
mation is in fact again the RPA form. 

For the self-energy, this approximation yields the form 

2(12) = /G(13)W(31), (3.27) 

the so-called GW approximation as introduced by He- 
din (1965). Of course, the Dyson equation for the 
Green's function G=G +G £G still has to be added, 
and even at the GW level one has to deal with a many- 
body self-consistent problem. However, the GW ap- 
proximation is a comparatively simple expression for the 
self-energy operator, which in principle allows the 
Green's function of an interacting many-electron system 
to be computed by starting from the Green's function 
G of a hypothetical independent-particle system with 
an effective one-electron potential, and iterate to self- 
consistency. The GW extends the well-known Hartree- 
Fock approximation, in which the self-energy is the ex- 
change potential X v , by replacing the bare Coulomb 
potential v by the dynamically screened potential W, 
e.g., X x = iGv is replaced by ~Z GW =iGW. The self- 
energy operator now consists of a dynamically screened 
exchange potential plus a dynamical Coulomb hole. In 
fact, the static approximation of the GW self-energy is 
none other than the COHSEX formula [Eq. (1.2)] dis- 
cussed in the introduction. GW has been shown to be 
physically well motivated and to be superior to the 
Hartree-Fock approximation, especially for metals, 
where Hartree-Fock (i.e., using the bare Coulomb po- 
tential) leads to unphysical results (Hedin and Lund- 
qvist, 1969). The dynamically screened interaction W in- 
troduces energy-dependent correlation effects absent in 
the one-particle picture. Further, an energy-dependent 
correlation decreases the Hartree-Fock band gap by 
raising the valence-band energy and lowering the 
conduction-band energy. There is some empirical evi- 
dence that supports the idea that even in the first itera- 
tion (that is, using just the noninteracting Green's func- 
tion G ) one obtains quite accurate results for one- 



electron properties such as the excitation energy 
(Hybertsen and Louie, 1985, 1986; Godby et al, 1986, 
1987, 1988; Aryasetiawan and Gunnarsson, 1998) and 
the quasiparticle lifetime (Campillo et al, 1999; Schone 
et al., 1999; Campillo, Pitarke, etal., 2000; Echenique 
et al., 2000; Campillo, Rubio, etal., 2000; Campillo, 
Silkin, et al, 2000; Keyling et al, 2000; Silkin et al, 2001; 
Spataru etal., 2001). This is important for practical ap- 
plications of the GW approach since, despite its formal 
simplicity, the practical solution of the self-consistent 
GW equations is a formidable task, which has been car- 
ried out only recently: self-consistent calculations were 
performed for the homogeneous electron gas (Holm and 
von Barth, 1998; Holm and Aryasetiawan, 2000; Garcia- 
Gonzalez and Godby, 2001), simple semiconductors, and 
metals (Shirley, 1996; Schone and Eguiluz, 1998). Self- 
consistency modifies the one-electron excitation spec- 
trum (excitation energies and lifetimes) as well as the 
calculated screening properties. The results turn out to 
be worse than those of the nonself-consistent G W 
calculations 18 or those obtained within the so-called 
GW , in which only the explicit G, and not W, is up- 
dated at every iteration, thus achieving partial self- 
consistency in G (von Barth and Holm, 1996). Improve- 
ments have been obtained, in particular concerning the 
position of plasmon satellites, using a partially self- 
consistent G, i.e. wave functions of order zero but up- 
dated quasiparticle energies (Bechstedt et al, 1994). The 
problem of self-consistency in G W calculations has been 
investigated more deeply in such simple systems as the 
homogeneous electron gas or an exactly solvable Hub- 
bard model. The main outcome of the self-consistent 
GW calculation for the homogeneous electron gas 
(Holm and von Barth, 1998; Holm and Aryasetiawan, 
2000; Garcia-Gonzalez and Godby, 2001) is that the total 
energy computed with the Galitskii and Migdal (1958) 
formula turns out to be strikingly close to the total en- 
ergy calculated using quantum Monte Carlo (Ceperley 
and Alder, 1980; Ortiz etal, 1999). Here few sum rules 
already determine most of the energy contributions in 
the homogeneous electron gas. Encouraging results are 
also obtained for the electron gas in three and two di- 
mensions, even for those ranges of densities for which 
the GW approach is often supposed to fail (Garcia- 
Gonzalez and Godby, 2001) and for inhomogeneous sys- 
tems such as bulk, surfaces, and interfaces (Sanchez- 
Friera and Godby, 2000; Garcia-Gonzalez and Godby, 
2002), and dimers (Aryasetiawan et al, 2002b). This re- 
sult may be related to the fact that the self-consistent 
GW scheme conserves electron number, energy, and to- 
tal momentum (that is, it fulfills the microscopic conser- 
vation laws; Martin and Schwinger, 1959, Baym and 
Kadanoff, 1961; Baym, 1962). However, concerning 
spectroscopic properties a systematic deterioration in 
the description of the bandwidth, quasiparticle excita- 



In this scheme G W , the full G has been replaced by G 
and the corresponding screened Coulomb potential is called 
W . 
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tion and lifetimes, and plasmon satellites is found. Re- 
sults along the same lines were obtained in the case of 
an analytically solvable finite Hubbard cluster recently 
used to compare the performances of several types of 
GW calculations (Verdozzi et al, 1995; Pollehn et al. 
1998; Schindlmayr et al, 1998). 

An important conclusion of this subsection is that, if 
one is interested in the excitation spectrum, self- 
consistent GW performs worse than the simpler G W 
scheme; efforts should be directed towards consistently 
including vertex corrections in the calculations (see next 
section). 

2. Second iteration step: vertex correction in P and 2 

Having arrived at the GW expression for the self- 
energy (or any approximation for X other than X = 0), 
one can go again through Eq. (3.22)-(3.25), starting with 
the latter. Now, for a nonvanishing functional derivative 
8%/8G, one obtains a correction to the bare vertex T 
= 8; this is the linear response of the self-energy to a 
change in the total potential of the system. This correc- 
tion affects the polarizability in Eq. (3.24) — the two 
Green's functions in P are no longer decoupled. In fact, 
vertex corrections account for exchange-correlation ef- 
fects between an electron and the other electrons in the 
screening density cloud. In particular this includes the 
electron-hole attraction (excitonic effects) in the dielec- 
tric response, which we shall look at later. 

The improved T and P can then be used, through Eqs. 
(3.23) and (3.22), in order to construct a new self-energy. 
This iteration step is still important, since, for example, 
long-range vertex corrections are needed in order to im- 
prove the structure of plasmon satellites (Aryasetiawan 
and Gunnarsson, 1998). Also the fact that self- 
consistency causes the results of GW spectroscopy cal- 
culations to deteriorate suggests that vertex corrections 
are important (Verdozzi et al, 1995). On the other hand, 
as noted by DuBois (1959), vertex corrections enhance 
correlation functions as the density-density response, 
whereas the inclusion of self-energy effects in the 
Green's function leads to a reduction. Similar conclu- 
sions were reached by Hong and Mahan (1994). There- 
fore one faces competing effects between self- 
consistency and vertex corrections: contributions from 
vertex corrections and self-consistency tend to cancel to 
a large extent for the 3D homogeneous electron gas. 
However, vertex corrections in quasi-2D electron-gas 
systems (like narrow quantum wells) are more impor- 
tant for the description of the quasiparticle band-gap 
renormalization and lifetimes (Marmorkos and Das 
Sarma, 1991). 19 The cancellation of vertex corrections 
with self-consistency seems to be a quite general feature. 
Vertex correction diagrams have also been found by de 



For the case of the semimetallic graphite (quasi-2D system) 
it has been shown by Spataru et al. (2001) that G W describes 
the observed energy dependence of the electron lifetimes but 
not the absolute value. 



Groot et al. (1996) to partially cancel the self- 
consistency effects in the case of a quasi-one- 
dimensional semiconducting wire. Their results show 
that, after systematic inclusion of all lowest-order cor- 
rections both to the RPA polarizability and to the GW 
self-energy, the band gap is roughly the same as at the 
first iteration G W . Cancellation of the first-order ver- 
tex and self-consistency corrections, however, has been 
found only to a minor extent for the band gaps of bulk 
silicon and diamond (Ummels et al, 1998). The authors 
have included vertex and self-consistency corrections to 
first order and found that corrections to the polarizabil- 
ity largely compensate each other, while this is not true 
for the gaps which become larger than the GW ones by 
0.7 and 0.4 eV for diamond and silicon, respectively. 20 

There are also cancellations between the same vertex 
correction used in different places; in fact, Mahan and 
Sernelius (1989) emphasize that a consistent way of han- 
dling vertex corrections is to include identical vertex cor- 
rections in both the polarization and the self-energy op- 
erator. 

The simplest improvement to the GW approximation 
in Eq. (3.22) consists in introducing a vertex correction 
consistent with the starting LDA calculation of the one- 
electron orbitals. The DFT exchange-correlation poten- 
tial is regarded as an approximation to the self-energy 
(Del Sole et al, 1994). Based on this idea, the vertex Y 
can be easily expressed in terms of the functional deriva- 
tive 8V xc ISp. This so-called GWY approximation 
(Rice, 1965; Hybertsen and Louie, 1986; Mahan and Ser- 
nelius, 1989; Del Sole et al, 1994; Mahan, 1994) reads in 
compact notation 



1 



1-fxcPo 



w=w l -= 

and thus 



q-fxcPo)v 

l-(v+f xe )P 



G=G 



X = iG W 1 T = iG 



l-(v+f xc )P - 



(3.28) 



(3.29) 



X is hence again of the GW form, but with an effective 
screening W 2 : =v/[l — (v+f xc )P ] which is in fact the 
testcharge electron dielectric function of linear-response 
theory (Del Sole et al, 1994; Hedin, 1999). The correc- 
tions given by f xc hence account for exchange- 
correlation effects at two levels. The first level concerns 
the electrons of the "screening medium," with an 
exchange-correlation hole around the electrons; when 
an electron is participating in the dielectric screening 



20 Note that there are also cancellations in the T-matrix ap- 
proximation where short-range interactions between two local- 
ized holes (electrons) can be described by a series of ladder 
diagrams exactly representing all repeated binary collisions 
(Kanamori, 1963). Due to the effects of the vertex corrections, 
the self-consistent version of the T matrix (using dressed 
propagators in the ladder) is found to be inferior to the non- 
self-consistent one (Cini and Verdozzi, 1986). 
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others are less likely to be found nearby. Second, the 
potential induced by the additional electron or hole also 
includes exchange-correlation interactions between the 
particle and the system electrons, and not only the clas- 
sical Coulomb interaction (since the particles are in fact 
indistinguishable). The vertex correction has to be put in 
both P and 2 in order to obtain this result (a vertex 
correction in P alone yields the testcharge-testcharge di- 
electric function). Results for silicon (Del Sole etal., 
1994) show that the quasiparticle energy gaps obtained 
using the RPA screening (as is done in standard GW 
calculations) or W 2 are close, whereas using 2 
= /G W 1 , i.e., putting the vertex correction only in P, 
alters the results, which gives evidence for the impor- 
tance of consistent corrections. Support for this sort of 
GWF approximation can already be found in the work 
of Rice (1965). He showed that, using Hubbard's energy 
expression corrected for exchange effects in an electron 
gas, a functional derivative of the total energy with re- 
spect to the electron occupation number gave a quasi- 
particle energy corresponding to t = iGoW 2 ■ 

The main drawback of the LDA-based GWF ap- 
proximation is the fact that results for bandgaps in semi- 
conductors as well as valence bandwidths are very close 
to the standard GW values. Further improvements 
might be obtained by using better approximations to the 
TDDFT kernel, including nonlocality and memory ef- 
fects (see Sec. V.C), or it might be necessary to work 
with the true nonlocality of the exchange-correlation 
self-energy, which means using a four-point kernel (see 
Sec. II. C). An example of such a calculation has been 
given by Shirley and Martin (1993), who performed GW 
calculations for atoms starting from Hartree-Fock and 
using a vertex correction consistent with the Fock ex- 
change term. Shirley and Martin found that the results 
for Group I, II, IV, and VIII elements are similar to 
those of ordinary GW without vertex corrections, again 
suggesting large cancellations. It should, however, be 
noted that the choice of the "reference" state was cru- 
cial. For example, electron removal energies from the 
neutral state were calculated as electron addition ener- 
gies to the corresponding singly ionized state. This might 
be seen as a partial inclusion of important vertex and 
self-consistency corrections already in the reference 
state, for calculations both with and without explicit ver- 
tex corrections. 

Systematic vertex corrections can hence be obtained 
through an iterative solution of Hedin's equations. One 
can push this scheme further by using the GW approxi- 
mation for 2 in the equation for the vertex function, and 
so on. An explicit formula for the vertex function corre- 
sponding to the full second cycle of iteration can be ob- 
tained which mixes certain diagrams of different orders 
in the screened interaction. Schindlmayr and Godby 
(1998) have shown that the vertex to order (n + l), 



can be expressed in terms of the lowest-order equiva- 
lents of all quantities except 2 . 



52 (n) (12) 
5G (,,) (45) 



xG ( " ) (46)G ( " ) (75)r(" + 1) (673), 



r (n + 1) (123)= 5(12)5(13) + J d(4567) 



r(" +1 )(123)= 8(12) <5(13) + J d(4567) 



<52 (,,) (12) 
5G (U) (45) 



x G (0) (46)G (0) (75)r (1) (673). (3.31) 

The appearance of the lowest order reduces the numeri- 
cal effort substantially. In particular G^ can be chosen 
such that it contains no satellite spectrum, but only a set 
of robust quasiparticle excitations. In the self-consistent 
limit the previous equation implies a relation be- 

tween the exact self-energy and vertex function. Nu- 
merical results for a model of strongly correlated elec- 
trons have indicated that this method yields improved 
excitation energies, in particular concerning the satellite 
spectrum. Nevertheless, to our knowledge this approach 
has never been applied to real systems. Instead, a 
slightly different version of the first iteration step is 
widely used, namely, the expression 

< + n f 52 (,,) (12) 

r(»+D(123) = 5(12) 5(13) + J ^(4567) g G („-i )(45) 



x G ( " ) (46)G ( " ) (75)r(" + 1) (673) 



(3.32) 



in its lowest-order version, 



r< 2 >(123) = 5(12)5(13)+ J~d(4567) 



52 (1) (12) 



(3.30) 



5G (U) (45) 

x G (1) (46)G (1) (75)r (2) (673). (3.33) 

Note that the only difference between Eqs. (3.30) and 
(3.32) is the version of the Green's function with respect 
to which the functional derivative of 2 is taken. The 
consequences might, however, be important, because 
whereas it is obvious that already the first-iteration so- 
lution of the integral equation (3.33) will shift the poles 
of P= — iGGT , this is not obvious in the case of Eq. 
(3.31). Rather, for n — 1 Eq. (3.31) multiplied by 
— iG^G^ yields a first-order polarizability with poles 
coming from the poles of both G^ and 

In fact, Eq. (3.33) is the starting point for the deriva- 
tion of the Bethe-Salpeter equation, which correctly 
yields features like bound excitons in the absorption 
spectra, as we shall show in Sec. IV.B.3. 

Obviously, choices like the one between Eqs. (3.30) 
and (3.32) are crucial only because one typically does 
not want to iterate to infinity, but to find the "best" first- 
iteration step, which is not well defined. This kind of 
problem regularly shows up when one is dealing with 
many-body Coulomb equations. 

Another way of partially summing higher-order dia- 
grams for describing vertex corrections to the valence 
electron Green's function is given by an exponential ex- 
pression very similar to the (almost exact) solution ob- 
tained for the core electron Green's function, as first 
shown by Hedin (1980). This approximation can also be 
derived as the first term in a cumulant expansion (Ar- 
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yasetiawan and Gunnarsson, 1998). The spectral func- 
tion in the cumulant expansion for a state below the 
Fermi level is expressed as 

A(k,co)=-^\ dte iMt e- ie k' +c ( k ''\ (3.34) 
2tt j _ m 

where n k is the occupation number of state k, e k is the 
single-particle eigenvalue, and C(k,t) is the cumulant 
operator, which can be derived in different ways (Hedin, 
1980). 21 The cumulant expansion contains only boson- 
type diagrams describing emission and reabsorption of 
plasmons; it does not contain diagrams corresponding to 
the interaction between a hole and particle -hole pairs. 
For this reason, the cumulant expansion primarily cor- 
rects the satellite description, whereas the quasiparticle 
energies are to a large extent determined by GW. This 
approximation has been applied successfully to describ- 
ing the multiple plasmon satellite structure in the spec- 
tra of alkali metals (Aryasetiawan et al. , 1996) as well as 
to substantially improving the calculated electron mo- 
mentum spectroscopy of graphite in regions of energy 
where no intensity is predicted from the LDA band 
structure (Vos et al, 2001). 

IV. GREEN'S FUNCTIONS IN PRACTICE 

From the previous subsection it emerges that the way 
quasiparticle energies and electron-hole excitations 
should be calculated is not unique. In practice schemes 
have been designed which allow one to obtain satisfac- 
tory results in many applications. An important variable 
in this strategy is the starting point of the calculations, 
i.e., one has to find the charge density (and often the 
lattice parameters and other ground-state properties) 
and some suitable one-electron eigenvalues and eigen- 
functions which allow construction of a starting one- 
particle Green's function. This task can be successfully 
carried out using the Hartree-Fock approximation in 
cases like atoms and molecules where screening is weak, 
or, most often, starting from static density-functional 
theory. In the following we shall summarize how the cal- 
culations are done in practice. 

A. Calculations of one-particle excited states 

First-principles GW calculations are today standard 
practice in many solid-state theory groups. Several re- 
views on quasiparticle calculations and the derivation 
and analysis of Hedin 's equations have been published 
in recent years (e.g., Aryasetiawan and Gunnarsson, 
1998; Aulbur et al. 1999; Hedin, 1999) and we refer the 



The simplest way is to identify the cumulant expan- 
sion {G(k,t) = ie- iBk ' +c(k -<'> = G (k,t)[l + C(k,t)+ 1 2C 2 (k,t) 
+ •••]} with the self-energy expansion of the Green's function 
(G = G +G XG +G 2G 2Go+ ••• )• In practice, using the 
GW approximation for the self-energy and a G determined 
from a DFT calculation, the cumulant is obtained equating 
G C = G 1G with t = l, GW -V xc . 



reader to those for further details. The fundamental re- 
view by Hedin and Lundqvist (1969) focuses on the elec- 
tron gas and contains a comprehensive introduction to 
all concepts and physical quantities and the connection 
to measurable spectra. In the recent reviews by Aryase- 
tiawan and Gunnarsson (1998) and Aulbur et al. (1999) 
the successes and shortcomings of the GW approxima- 
tion and numerical implementations (plane waves, local 
orbital basis, real-space/imaginary-time) and tricks to 
improve the computational efficiency are described, to- 
gether with a detailed compilation of results obtained 
for semiconductors, transition-metal oxides, fullerenes, 
superlattices, interfaces, surfaces, defects, Schottky bar- 
riers, and atoms and molecules. Apart from their impor- 
tance on their own, quasiparticle energies also serve as 
an input to further spectroscopy calculations, in particu- 
lar for the Bethe-Salpeter equation approach. Therefore 
in the following we outline the way GW calculations are 
usually performed and give some illustrative examples. 

1. GW calculations 

A GW calculation requires in principle solution of 
Eq. (3.16) for the quasiparticle energies and amplitudes, 
taking for X the product of the Green's function G and 
the screened Coulomb interaction W. The complexity of 
GW leads, in practical applications, to approximations 
in the construction of G and W, which we shall examine 
in the following. 

a. Evaluation of the Green's function G 

In practice G is constructed using single-particle or- 
bitals from a Kohn-Sham (or, less often, from a Hartree- 
Fock) calculation. One considers the single-particle or- 
bitals and eigenvalues as a zeroth-order approximation 
to the quasiparticle amplitudes and energies (see Hedin, 
1995, for a discussion of this issue. This suggests the pos- 
sibility of looking for a first-order, perturbative solution 
of Eq. (3.16) with respect to (X — V xc ) or to (X — X x ), 
where V xc is the exchange-correlation potential of 
Kohn-Sham (or X x is the exchange potential of Hartree- 
Fock). This is the approach followed in many practical 
GW calculations for real systems. 22 The diagonal matrix 
elements of (X — V xc ) give the quasiparticle energies as 

ef -£j +Z i {(f> j |X(e,- )-V xc \(p i ), (4.1) 

with Z~r 1 = l-{<f>f s \d%/de\ e Ks\</>f s ). In the case of 
simple bulk semiconductors, this approximation turns 
out to be very close to the exact result obtained by di- 
agonalizing H H +% GW (but without any self-consistent 
update of G and W; Hybertsen and Louie, 1986). In 
other cases, such as transition metals and finite and low- 
dimensional systems, the validity of the perturbative ap- 
proach in X — V xc should be checked explicitly. For in- 



The solution of the exact local Kohn-Sham exchange has 
also been used as input for the GW quasiparticle calculation 
(Aulbur et al, 2000). The quasiparticle corrections seems to be 
smaller in this case. 
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stance, Pulci et al. (1999) have shown that for the (110) 
surface of GaAs, even if quasiparticle energies are not 
altered when doing the full diagonalization of the quasi- 
particle Hamiltonian, the quasiparticle wave functions 
corresponding to nearly degenerate states, with different 
degrees of surface/bulk localization, strongly hybridize, 
making the quasiparticle wave functions significantly dif- 
ferent from the LDA ones. Differences between LDA 
and GW wave functions have also been found in the 
case of SiH 4 for states crossing the level of zero poten- 
tial (Rohlfing and Louie, 2000). In general, differences 
between LDA and quasiparticle wave functions can at 
least be expected in any nonbulk system, i.e., when there 
are regions of space where the density goes to zero. In 
these regions, LDA produces an exchange-correlation 
potential with an incorrect asymptotic behavior (e.g., 
missing the — II r tail in the case of atoms and clusters); 
hence, LDA wave functions may have a wrong spatial 
localization. The GW potential, on the other hand, does 
have the correct asymptotic behavior (Charlesworth, 
et al, 1993; Garcia-Gonzalez and Godby, 2002). 23 

b. Evaluation of the screened Coulomb interaction W=e~ V 

One of the main difficulties is the evaluation of the 
full dielectric response of the system. In fact, even within 
the chosen zeroth-order scheme, i.e., within RPA-LDA, 
the calculation of e _1 (r,r',a)) remains a difficult task 
from a numerical point of view (due to the spatial non- 
locality and frequency dependence of e). For this rea- 
son, approximated schemes have been developed: 

• The model dielectric function (Levine and Louie, 1982; 
Hybertsen and Louie, 1988; Bechstedt, Del Sole, 
et al, 1992). These approaches allow one to consider- 
ably reduce the computational effort or even to avoid 
completely the ab initio calculation of e. 

• The plasmon-pole approximation (Hybertsen and 
Louie, 1986; von der Linden and Horsch, 1988; Ha- 
mada et al., 1990; Engel and Farid, 1993). Most mod- 
els are based on the observation that e(G,G',w) _1 is 
generally a peaked function in co that can be approxi- 
mated by a single-pole function in co. The pole posi- 
tion and strength are determined by imposing sum 
rules (Hybertsen and Louie, 1986), or by fitting each 
element e _1 (G,G') at two points co along the imagi- 
nary energy axis (Godby and Needs, 1989). Since the 
evaluation of £ involves an integration over the en- 
ergy, the fine details of the energy dependence are not 
critical, and the plasmon-pole approximation turns 
out to work reasonably well. One important drawback 
is that quasiparticle lifetimes cannot be computed 
within a plasmon-pole approximation, as Im(X) is 
zero everywhere except at the pole. The validity of 
the plasmon-pole approximation has not been system- 



atically tested, except in a very few cases (see, for 
example, Aulbur et al., 1999; Arnaud and Alouani, 
2000). 

When the full energy dependence of the dielectric ma- 
trix is retained, the integration in co may be performed 
along the imaginary axis, where e _1 is well behaved 
(Godby et al. 1988; Schone and Eguiluz, 1998), by pick- 
ing up all the poles along the real axis (Aryasetiawan, 
1992; Fleszar and Hanke, 2000), or by using the 
transition-space spectral representation of s _1 , which 
allows one to perform the frequency integration analyti- 
cally (Shirley and Martin, 1993). 24 



c. Application of the self-energy 2 

For many applications, the simple prescriptions de- 
scribed above have yielded results within 10-15 % of the 
experimental ones, the typical example being the direct 
quasiparticle gap of diamond, as discussed in the intro- 
duction. An agreement of the same quality between 
quasiparticle energies and photoemission or inverse 
photoemission data has also been obtained for more 
complex systems like surfaces, atoms, and nanostruc- 
tures (Aulbur et al, 1999), and the GW calculations are 
systematically used to study quasiparticle excitations in 
realistic systems of practical interest. For example, in the 
case of an ethylene molecule (C 2 H 4 ) adsorbed on the 
Si(001)-(2x 1) surface GW seems to improve the cal- 
culated tunneling currents measured in scanning tunnel- 
ing microscopy (Rignanese et al, 2001). 25 A recent GW 
study on YH 3 illustrates the important consequences of 
self-energy corrections (Miyake et al. , 2000; van 
Gelderen et al, 2000). It turns out that GW corrections 
remove the band overlap responsible for erroneous me- 
tallic LDA behavior in YH 3 that has a measured optical 
gap of 2.8 eV. Hence the gap is of electronic origin 
rather than structural, and not due to strong correlation 
(Eder et al, 1997; Ng et al, 1997). 

Another physical quantity that can be obtained from a 
knowledge of the full, complex self-energy 2 is the qua- 
siparticle lifetime, i.e., the electron-electron scattering 
contribution to the linewidth. The damping rate of an 
excited electron in the state ^o( r ) with energy E is in 
fact obtained as (Echenique et al. , 2000) 



23 One recent G W application touching this point is the work 
of White et al. (1998) on the image potential at the Al(lll) 
surface. 



24 Mixed-space (Blase et al, 1995) and real-space/imaginary- 
time techniques (Rojas et al, 1995), as well as a technique to 
eliminate the unoccupied state summations (Reining et al. , 
1997), can also be used to reduce the computational effort (see 
the reviews by Aulbur et al. , 1999; and Aryasetiawan and Gun- 
narsson, 1998, for a detailed comparison of different numerical 
implementations) . 

2 Similar studies concerning the interpretation of scanning 
tunneling microscopy images have also been performed in the 
framework of the LDA+U approach [see, for example, the 
work by Dudarev et al. (1997) on the NiO(100) surface]. We 
refer the reader to Anisimov et al. (1997) for a description of 
the relation between GW and LDA+U. 
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FIG. 2. Comparison of the calculated and experimental band 
structures for copper: solid line, GW quasiparticle excitation 
energies (Marini et al, 2002); dashed line, DFT-LDA eigenval- 
ues; O, experimental data compiled by Courths and Hiifner 
(1984). A comparison to more recent experimental data (Stro- 
cov et at, 1998, 2001) yields the same agreement. 



t- 1 =-2 J dr j dr' ifj*(r)lmX(r,r\E)Mr'). (4.2) 

In the GW approximation, the imaginary part of 2 is 
determined by the imaginary part of e _1 that contains 
all decay processes of the initial state. Quasiparticle 
damping times were calculated for diamond by Strinati 
et al. (1980, 1982), who found results consistent with 
photoemission data. A treatment of the electron dynam- 
ics, including band structure and dynamical screening ef- 
fects, is necessary for quantitative comparisons with ex- 
periment (see, for instance, Biirgi et al, 1999; Valla et al, 
1999b). An illustrative example is given by Campillo 
et al. (1999; Campillo, Pitarke, et al, 2000; Campillo, Ru- 
bio, et al, 2000) and Gerlach et al. (2001), who evaluated 
the quasiparticle lifetimes of electron and holes in bulk 
Cu within the GW scheme, showing an increase in the 
lifetime close to the Fermi level as compared to the pre- 
dictions of a free-electron-gas model of the solid. How- 
ever, the lifetimes get closer to those of a free-electron 
gas for hole energies below ~ 3 e V, whereas the experi- 
mental data show a distinct asymptotic behavior. This 
discrepancy between theory and experiment may be a 
signature of important departures of the DFT-LDA 
band structure of Cu, which is used in the calculation of 
W, from the actual quasiparticle band structure (see Fig. 
2; Marini et al, 2002). 

GW results for Cu agree with photoemission data 
within 30 meV for the highest d band, correcting 90% of 
the LDA error. The energies of the other d bands 
throughout the Brillouin zone are reproduced within 300 
meV, and the maximum error ( — 600 meV) is found for 
the bottom valence band at the T point, where only 50% 
of the LDA error is corrected. This level of agreement 
for the d bands cannot be obtained without including 
self-energy exchange contributions coming from 3 s and 



3p atomic levels, 26 demonstrating the importance of 
those contributions. On the other hand, the total band- 
width is still larger than the measured one. This overes- 
timate of the GW bandwidth with respect to the experi- 
mental one seems to be a quite general feature, which is 
not yet properly understood. (See Yasuhara et al, 1999, 
for a discussion about Na, Ku et al , 2000 for a comment, 
and Yasuhara et al, 2000 for reply.) 

Fully self-consistent GW calculations have been per- 
formed by Schone and Eguiluz (1998) for crystalline sili- 
con and potassium, with the conclusion that for real ma- 
terials, self-consistent GW without vertex corrections in 
X yields quasiparticle energies and bandwidths that dis- 
agree with experiment; for example, the absolute gap of 
Si turns out to be 1.91 eV, compared with the experi- 
mental 1.17 eV. On the other hand, Massidda et al. 
(1995) have calculated self-consistently the electronic 
structure of MnO, using a frequency-independent model 
£ derived by Gygi and Baldereschi (1989). An encour- 
aging agreement was found with experiment concerning 
the energy gap, magnetic moment, bandwidth, and spec- 
tral distribution of Mn 3d states. The fact that self- 
consistency leads to poorer results in the calculations of 
Schone and Eguiluz (1998), but leads to an improvement 
in the calculations of Massidda et al. (1995) shows that 
the dynamical aspect of the self-energy is very problem- 
atic in this context. 

Few works including vertex corrections exist to date 
on real materials. One example is given by the GWT 
calculation on silicon, described in Sec. III.C.2, where T 
is taken from a DFT-LDA approach. However, the ef- 
fects of this approximate kernel are rather small. It is 
difficult to guess what the effect would be if a better (in 
particular more long-range) approximation of the kernel 
were used. The work of Ummels et al (1998), discussed 
in the same section, suggests that some changes might be 
found. 

The effects of vertex corrections can also explain why 
GW calculations have not always been successful in de- 
scribing self-energy effects, for example for ^-electron 
bands (where the problem might come from strong cor- 
relation or simply from self- interaction effects, i.e., from 
the strong localization). A critical discussion about the 
extent to which the GW approximation is capable of 
describing highly correlated systems such as NiO is 
given by Aryasetiawan and Gunnarsson (1995). 

Computationally, the evaluation of 2 as G W is a hard 
task. 27 Calculations can scale as badly as N 4 , where N is 
the number of atoms, hence reaching the limit of com- 
puter power well before DFT-LDA calculations, which 
usually scale as TV 3 . 



26 The relevance of semicore states in the self-energy was also 
pointed out by Rohlfing et al. (1995, 1998). 

27 A typical bottleneck is given by the summations over the 
empty states (Reining et al, 1997), both for the determination 
of the screening and for G . Another difficulty comes from the 
convergence of Brillouin zone integrals (i.e., k -point sampling) 
of functions that have a factor |k+G| 2 in the denominator as 
the exchange term (Pulci, 1998). 
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It turns out that, whereas the "standard" prescription 
for GW calculations yields good results for many appli- 
cations, one has to be careful when effects beyond this 
perturbative-RPA scheme have to be included. In order 
to gain some insight concerning the choices to be made, 
it is useful to resort to the outcome of model calcula- 
tions. 

2. Comparison of GW and DFT results 

As mentioned in Sec. II. A, the Kohn-Sham eigenval- 
ues, appearing as Lagrange multipliers in the minimiza- 
tion of the density functional, do not carry a direct 
physical meaning, even though they can be considered 
as well-defined approximations to the excitation ener- 
gies (Gorling, 1996; Filippi, 1997). In particular, they can 
be viewed as a good zero-order starting point for pertur- 
bative self-energy calculations (see Sec. IV.A.l). In fact, 
there is no equivalent of Koopman's theorem for the 
Kohn-Sham eigenvalues, and a direct comparison of the 
LDA energy gap of a semiconductor or insulator with 
the experimental (photoemission and inverse photo- 
emission) band gap normally yields a severe underesti- 
mate, of the order of 30-50 % or more (see, for ex- 
ample, Bechstedt, 1992). This error is corrected in a GW 
calculation of quasiparticle energies. 

For a more precise and general comparison between 
DFT and GW, one has to distinguish between the ef- 
fects of the approximated exchange-correlation func- 
tional, for example LDA (which affects the total energy 
differences and the effective potential, hence the eigen- 
values) and the effect of interpreting Kohn-Sham eigen- 
values as electron addition and removal energies (see, 
for example, Jones and Gunnarsson, 1989). Almbladh 
and von Barth (1985) showed that in exact DFT the 
highest occupied eigenvalue (highest occupied molecu- 
lar orbital level, or HOMO) does have a physical inter- 
pretation: it corresponds to the ionization potential. 
Hence, in exact DFT, the ionization potential could be 
derived directly from a single calculation for the ground 
state of an iV-electron system. Unfortunately, the use of 
approximated exchange-correlation functionals (like 
LDA) quite severely affects the Kohn-Sham highest oc- 
cupied level. Hence, when comparing excitation ener- 
gies directly with Kohn-Sham eigenvalues, one must dis- 
tinguish two cases: the highest occupied level, for which 
the error is due only to the LDA, and the other levels 
[e.g., the lowest unoccupied molecular orbital (LUMO) 
eigenvalue], for which, even in exact DFT, a discrepancy 
could be found. For example, in the Be atom (Jones and 
Gunnarsson, 1989), both errors are large: the LDA 
HOMO ( — 5.6 eV) misses the ionization potential 
(-9.32 eV) by almost 4 eV, and the HOMO-LUMO gap 
is 3.5 eV in LDA, 3.62 eV in exact DFT, and larger than 
9 eV experimentally. Part of the error comes from the 
spurious self-interaction which is contained in the LDA. 
Perdew and Zunger (1981) have proposed a self- 
interaction-corrected approach in order to eliminate this 
contribution. This is also the case in realizations of DFT 
which allow the inclusion of the exact local exchange 



potential as an optimized effective potential (OEP) (de- 
fined as the functional derivative of the exact exchange 
energy E x with respect to the electron density) and 
which have been introduced with great success (Gross 
et al., 1996; Gorling et al, 1999). 

On the other hand, in a finite system it is possible to 
compute quasiparticle energies directly as energy differ- 
ences, provided one is able to compute the total energy 
of the system with N electrons in the ground state, and 
that of the system with N±l electrons in the ground 
state or some excited state. This scheme, known as the 
delta-self-consistent-field method, has been successfully 
used within DFT to describe bound one-electron excita- 
tions and the ionization potential of atoms and mol- 
ecules (Jones and Gunnarsson, 1989). 

Hence, DFT- or Hartree-Fock-based delta-self- 
consistent-field calculations can be directly compared 
with GW results. For atoms in the iron series, Shirley 
and Martin (1993) have shown that GW yields results of 
a quality similar to, or slightly worse than, those com- 
puted from Hartree-Fock or local-spin-density total en- 
ergy differences (Baroni, 1984; Vukajlovic et al., 1991). 
Small sodium clusters can be taken as another example: 
Na 4 and Na 6 vertical and adiabatic ionization potentials 
obtained with the delta-self-consistent-field approach 
(Martins et al., 1985) compare well with DFT+ GW re- 
sults for jellium spheres (Saito et al., 1989) and for the 
real clusters (Onida et al., 1995; Reining et al., 2000). 
Similarly, Ogiit et al. (1997) studied the evolution of the 
quasiparticle and optical gaps of hydrogenated silicon 
clusters as a function of cluster size. Delta-self- 
consistent-field LDA calculations were used to estimate 
the quasiparticle correction to the HOMO-LUMO gap. 
Ogiit et al. concluded that quantum confinement, as well 
as reduction of screening due to finite size, leads to ap- 
preciable excitonic corrections. This correction turns out 
to be large and size dependent. The delta-self-consistent 
field results point out that for finite systems the Hartree 
relaxation of the wave functions when the electron is 
added or removed is essential, while it is of no relevance 
for delocalized states in bulk solids, as pointed out by 
Godby and White (1998; see also the reply of Ogiit et al, 
1998). In infinite systems, the correlation part takes over 
in relevance (see the discussion in Sec. VI. C). 

Finally, another way to overcome the "gap problem" 
without performing a self-energy calculation has been 
proposed by Mackrodt et al. (1996): Hartree-Fock calcu- 
lations on nickel oxide (NiO) and on the lithium-doped 
material have allowed the determination of the band 
gap of the former, by taking eigenvalue differences be- 
tween two empty states of the latter. Simply taking the 
difference between Hartree-Fock eigenvalues of an oc- 
cupied and an empty state of NiO would lead to far too 
large a gap. 

In conclusion, Kohn-Sham eigenvalues are generally 
in disagreement with photoemission and inverse photo- 
emission energies, both for finite and for extended sys- 
tems, with errors of 30-50 % or more. On the other 
hand, self-energy-corrected quasiparticle energies (ob- 
tained via the GW method) usually agree well with pho- 
toemission and inverse photoemission data: the reported 
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errors are typically in the range of 10-15 % of the band 
energies measured with respect to the Fermi level, hence 
being in absolute value larger for systems with a larger 
gap. In the case of finite systems, the delta-self- 
consistent-field method also yields accurate excitation 
energies, since wave-function localization allows for fi- 
nite relaxation effects. In infinite systems with delocal- 
ized wave functions the delta-self-consistent-field 
method is inadequate, since Hartree relaxation effects 
vanish and correlation must be taken into account. GW 
contains both relaxation and correlation and is therefore 
suitable to describe electron addition and removal ener- 
gies. 

In the next section we shall discuss the use of Green's- 
function theory for the study of two-particle excited 
states, such as the electron-hole pairs that are created in 
absorption measurements. For an overview of earlier 
and some more recent work, see Hanke et al. (1983) and 
Rohlfing and Louie (2000), respectively. 

B. Calculations of two-particle excited states 

The key quantity in the calculation of two-particle ex- 
cited states is the polarizability P, as introduced in Eq. 
(1.3). 

1. Independent-quasiparticle approximation 

As discussed in Sec. III.B, P is given by Eq. (3.24). In 
Sec. III.C we defined the RPA form P=-iGG [Eq. 
(3.26)], from which the screened interaction in GW cal- 
culations is usually obtained. In principle, Eq. (3.26) 
could also be used for the calculation of absorption spec- 
tra and other spectra that involve the creation of 
electron-hole pairs. P is then just the product of two 
one-particle Green's functions. Using Eq. (3.4) and after 
a Fourier transformation from time to frequency space, 
one obtains for the independent-(quasi)particle P the 
result (1.4) in its time-ordered form, from which the re- 
tarded version can be deduced. As in the case of the 
one-particle Green's function, the irreducible polariz- 
ability P (which is derived from a two-particle Green's 
function) can be understood in its Lehmann representa- 
tion or in the complex plane. In the second case, which is 
the representation used in practical applications, the 
electron and hole energies entering the denominator can 
be complex, leading to a finite lifetime of the two- 
particle excitation (because of the finite lifetime of the 
electron and/or the hole). Additional deexcitation chan- 
nels are included when the electron-hole interaction is 
taken into account; this will be discussed below. 

Formula (1.4) is at present still the standard expres- 
sion used for many ab initio calculation of optical spec- 
tra of real materials. Note that, whereas the form of P is 
prescribed by Eq. (3.26), the iteration approach is again 
not uniquely defined concerning the ingredients. Hence 
the Green's functions entering Eq. (3.26) may be consid- 
ered to stem, for example, from a Kohn-Sham or from a 
GW calculation. In practice, most often DFT Green's 
functions are used when the approximation P = P I q P 
[see Eq. (1.4)] is made. This can be understood as an 



approximation to the Green's-function iteration scheme, 
but also as an approximation to TDDFT, as will be ex- 
plained in Sec. V. Sometimes a mixed DFT- G W Green's 
function is used; the eigenvalues entering G are taken to 
be GW ones, whereas the wave functions are Kohn- 
Sham orbitals. This approach (called GW-RPA in the 
introduction) was first tried in order to overcome the 
redshift of DFT-based absorption spectra with respect to 
experiment, but in general it overcorrects, and, more- 
over, does not improve line shapes. 28 Today, GW-RPA is 
mostly seen as a good approximation to the first- 
iteration result for P jq P , which is then used to include T 
in the second iteration, i.e., it is the starting point of the 
Bethe-Salpeter approach, which will be outlined in the 
following. We shall concentrate on absorption spectra, 

1. e., on the calculation of P [Eq. (2.9)], where going be- 
yond the RPA gives the most visible effects. All results 
can be easily generalized to the case in which one is 
instead interested in ^, replacing the local field effects 
contribution v by v (see Sec. II. B). 

2. Electron-hole attraction 

The inclusion of vertex corrections [i.e., the inclusion 
of r in Eq. (3.24) for P] can be achieved through a sec- 
ond iteration of Hedin's equations. This means that one 
has to go again from Eq. (3.22) to Eq. (3.25), now using 
~% = iGW in the latter. This yields an integral equation 
for r, 

r(123) = 8{12) <S(13) 

+ /W(l + 2) J d(67)G(16)G(72)r(673). 

(4.3) 

Here, the approximation S%l 8G = iW is used. This 
means that (i) one neglects the term iG(SW/SG), 
which contains information about the change in screen- 
ing due to the excitation and is considered to be small, 
and that (ii) Eq. (3.32), and not Eq. (3.30), has been 
used. As discussed above, this latter choice in the itera- 
tion scheme can be crucial when only one or a few itera- 
tions are performed. 

One can transform Eq. (4.3) to an integral equation 
for a generalized 3 P, defined as 

3 P(312)=-/J" d(67)G(16)G(72)r(673), (4.4) 

by multiplying with — «G(41)G(25) on the left and in- 
tegrating over d(l2): 

3 P(345)=-/G(43)G(35) 

+ ij <i(12)G(41)G(25)W(l + 2) 3 P(312). 

(4.5) 



It does, however, yield rather good static dielectric con- 
stants (Levine and Allan, 1989). For silicon and germanium, 
the error decreases by one order of magnitude when GW in- 
stead of LDA eigenvalues are used. 
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FIG. 3. Feynman diagrams representing the Bethe-Salpeter 
equation for x- 



The polarization P can hence in principle be obtained 
from P(31) = 3 P(311). However, the kernel GGW of 
the integral equation (4.5) is a four-point function. 
Therefore, as outlined in Sec. II. C, it is equally conve- 
nient to introduce a four-point 4 P, which has the advan- 
tage that local field effects can be included on the same 
footing, i.e., one can directly obtain P (or x). In fact, 
closely following the prescriptions in Sec. II. C, one can 
introduce the four-point dynamically screened interac- 
tion 



1 W( 1234) = W( 12) S( 13) S( 24) . 



(4.6) 



Then one gets for Eq. (4.5) the four-point integral equa- 
tion 



4 P = *P IQP -*P 4 W 4 P IQP , (4.7) 

with the obvious generalization from Pjq p to 4 P jq P [see 
Eqs. (1.4) and (B16)]. One can now also write Eq. (2.10) 
for four-point quantities, namely, 



4 P( 1234) = 4 P( 1234) 

+ J d(5678) 4 P( 1256) <5(56) S(78)v (57) 

x 4 P(7834). (4.8) 

Putting these equations together, one obtains the Bethe- 
Salpeter equation for A P: 



4 P = 4 P IQP + 4 P IQP K 4 P, 



(4.9) 



where the kernel K contains an electron-hole exchange 
contribution v and electron-hole attraction — W: 

K( 1234) = S( 12) <?(34)u ( 13) - S( 13) 5(24) W{ 12) . 

(4.10) 

The equation for the corresponding 4 x differs only by 
the fact that v instead of v appears. The diagrams rep- 
resenting the Bethe-Salpeter equation for x are shown 
in Fig. 3. The Green's-functions lines are dressed, i.e., 
they are quasiparticle lines. 

One can make an immediate connection to the discus- 
sion of the four-point equations in Sec. II. C. In fact, 



most often W is taken to be the statically screened Cou- 
lomb interaction. Equation (4.10) for K corresponds 
then to the time-dependent screened Hartree-Fock 
approximation. 29 In particular, one can follow the 
scheme of Sec. II. C and Appendix B.2, and derive the 
effective two-particle Hamiltonian H 2p defined in Eq. 
(2.22), which is often used in the context of the Bethe- 
Salpeter equation approach. In fact, when static screen- 
ing is used in W the effective two-particle equation takes 
the particularly simple, energy-independent form, 

-^^p^^mI" 1 " 2 ' (4-ii) 



where it is understood that matrix elements of the ker- 
nel (4.10) are taken with respect to four one-particle 
orbitals n 1 ---n 4 : U( nin2 )(„ 3 „ 4 ) = («in 2 |5'l"3«4> and 
^(n 1 « 2 )(n 3 n 4 ) = ("i w 3l W\n 2 n 4 ). The solution of this 
equation allows one to construct then the absorption 
spectrum from Eq. (2.23). If one considers only the reso- 
nant part of H 2p , i.e., the part mixing only transitions of 
positive frequency, the resulting operator is Hermitian 
and its eigenstates orthogonal. Therefore one obtains 
the simpler formula 



e M (a») = l-limi;o(q)2 



|2 ni „ 2 <ni|e- iq > 2 M 
co — E k + irj 



«1"2|2 



(4.12) 



As discussed below, this turns out to be a very good 
approximation for the calculation of bulk absorption 
spectra. 

Summarizing this section, we recall that (i) the Bethe- 
Salpeter approach to the calculation of two-particle ex- 
cited states is a straightforward extension of the GW 
approach for the calculation of one-particle excited 
states; and (ii) it leads to an effective two-particle 
Hamiltonian H 2p which is of the general form found in 
time-dependent problems, as explained in the introduc- 
tion. The electron-hole interaction that appears in H 2p 
has two contributions: the first involves the bare Cou- 
lomb interaction (v or v , depending on whether x or P 
are calculated), which connects the electron and hole 
indexes in an exchangelike manner and is therefore also 
called electron-hole exchange. In order to avoid confu- 
sion, however, it should not be forgotten that this con- 
tribution stems from the density variation of the Hartree 
term and contains just the local field effects in the case 
of v. The second contribution involves W, which con- 
nects the electron and hole indexes like a direct Cou- 



It corresponds only for the kernel K, because the eigenval- 
ues entering Pjqp are in practice calculated in the full (dy- 
namic) GW scheme, and the eigenfunctions are approximated 
with LDA ones, not screened Hartree-Fock ones. 
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21omb interaction and is therefore called the screened 
electron-hole interaction. It comes from the variation of 
the self-energy, in other words, from an exchangelike 
term in the one-quasiparticle Hamiltonian. It is this last 
term that is responsible for the appearance of bound 
states, and it can even lead to hydrogenlike spectral fea- 
tures in insulators. 

In practice, the spin a n is most often not considered 
explicitly. One then makes use of the fact that only (u,c) 
transitions with cr v = cr c (i.e., singlet excitons) contribute 
to the optical spectrum. Since W does not contain a spin 
interaction, and because of the way the S functions con- 
nect indexes, the matrix elements Wr vc \/ tt i c i\ are diago- 
nal in spin, i.e., a v = a v i . Instead, the matrix elements 
V( vc ){v'c') are independent of a v and cr v , . Hence spin is 
implicitly included via a factor of 2 for v in Eq. (4.11). 



3. The effects of the electron-hole interaction 

Exciton calculations for real systems within a fully 
first-principles scheme, starting from DFT, have recently 
been performed. However, the importance of the exci- 
tonic effects has been known for a long time, even in 
cases in which the features are less characteristic than, 
the appearance of a hydrogenic series of bound states in 
the band gap. Excitonic effects in the absorption line 
shape of semiconductors above the fundamental gap 
were calculated some time ago (Hanke and Sham, 1974, 
1975, 1980; del Castillo and Sham, 1985). Hanke and 
Sham (1980) solved the Bethe-Salpeter equation [Eq. 
(4.9)] for the particle-hole response function of bulk sili- 
con, using a linear combination of atomic orbitals basis, 
with a semiempirical band structure fitted to optical ex- 
periments. They found important corrections to the 
independent-particle result, in particular a strong in- 
crease in the lowest main absorption peak (El), leading 
for the first time to a calculated absorption spectrum of 
bulk silicon in qualitative agreement with experiment. 

The ab initio approaches used today for solution of 
the Bethe-Salpeter equation (Albrecht et al, 1998a, 
1998b; Benedict et al, 1998a; Rohlfing and Louie, 
1998b) mostly follow the scheme introduced by Onida 
et al. (1995) for the spectrum of the cluster Na 4 and that 
of Albrecht et al. (1997) for the optical gap of bulk 
Li 2 0. This procedure consists of (i) a ground-state DFT 
calculation; (ii) a GW calculation to correct the eigen- 
values; and (hi) solution of the Bethe-Salpeter equation 
using GW eigenvalues, Kohn-Sham orbitals, and static 
RPA screening for the electron-hole interaction. It thus 
includes several steps and choices, which are schemati- 
cally represented by the flow diagram in Fig. 4. 

This scheme, using the approach of the effective two- 
particle Hamiltonian, has again been applied to the cal- 
culation of the absorption spectrum of bulk silicon (Al- 
brecht et al., 1998a), showing that the Bethe-Salpeter 
method used without any adjustable parameter can yield 




FIG. 4. Flow diagram sketching the solution of the Bethe- 
Salpeter equation in practice. 

absorption spectra of continuum excitons in bulk semi- 
conductors in quantitative agreement with experiment 
(see Fig. 5). 30 

On the other hand, as pointed out above, the electron- 
hole attraction can also lead to bound excitons in insu- 



The present curve has been calculated by Olevano and 
Reining (2000) using an improved Brillouin-zone sampling 
with respect to the original publication of Albrecht et al. 
(1998a). See also the discussions by Cardona et al. (1999) and 
Albrecht et al. (1999). 
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FIG. 5. Silicon absorption spectrum [Im(£ M )]: •, experiment 
(Lautenschlager et al., 1987); dash-dotted curve, RPA, includ- 
ing local field effects; dotted curve, GW-RPA; solid curve, 
Bethe-Salpeter equation. 

lators, where the interaction is only weakly screened. 
The resulting spectra can also be well described by the 
ab initio Bethe-Salpeter approach, as has been illus- 
trated by the work of Benedict et al. (1998a). These au- 
thors use a recursive method (Haydock, 1980) to invert 
the Bethe-Salpeter equation (4.9), which has the advan- 
tage that one can make use of the <5 functions in the W 
contribution to the kernel in order to make calculations 
less cumbersome. 31 Figure 6 shows the result for bulk 
LiF (Shirley, 2001). 32 The sharp peak in the absorption 
spectrum at about 12 eV is clearly due to excitonic 
effects — the spectrum without excitons is essentially fea- 
tureless in the range of the absorption onset. The 12-eV 
peak defines the optical gap, which turns out to be only 
about 0.5 eV lower than the experimental one, whereas 
the independent-quasiparticle result shows a consider- 
able overestimate. 

The same approach has also been used by Benedict 
et al. (1998b) to obtain the absorption spectra of bulk Si, 
C, Ge, and GaAs. Moreover, they have studied the 
wide-gap semiconductor GaN (in its wurtzite and zinc- 
blende structures), and the insulating CaF 2 crystal 
(Benedict and Shirley, 1999). In all cases, the inclusion 
of excitonic effects turned out to be crucial for a quan- 
titative agreement between theory and experiment. 
Rohlfing and Louie (1998b) have introduced a different 
way of extending the applicability of the Bethe-Salpeter 



jl The calculations are simplified further by the use of a model 
electron-hole screening. 

l2 Shirley's calculation is the same as that of Benedict et al. 
(1998a), but with an improved Brillouin-zone sampling and 
including more bands. 
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FIG. 6. LiF, absorption spectrum [£2 = I m ( £ M)]- Continuous 
curve with crosses, experimental (Roessler and Walker, 1967). 
Dashed curve, Eq. (2.8) using GW eigenvalues. Plain continu- 
ous curve, Bethe-Salpeter equation (Shirley, 2001). 

approach: they noticed that one of the bottlenecks in the 
calculation of bulk spectra is the determination of the 
(N v XN c XN k ) 2 matrix elements of H 2 { l kck)(v , k , c , kl) , 
where N v and N c are the numbers of valence and con- 
duction bands that have to be taken into account in the 
construction of the spectrum, and N k is the number of k 
points at which the direct vk^ck transitions are consid- 
ered. They hence proposed an interpolation scheme in k 
for the H/ ', .. u , , , .,, M , which allowed them to calcu- 
late the absorption spectrum of bulk GaAs including the 
extremely weakly bound exciton states (some meV of 
binding energy) in the gap. The description of the latter 
feature in fact required the use of about 1000 k points in 
a sphere of 0.015-a.u. radius around k = 0. 

Rohlfing and Louie (1998b) and Benedict and Shirley 
(1999) also noticed that the excitonic effects above the 
gap are due to the mixing of transitions given by the 



coefficients A 



(«i«2) ; 



in Eq. (4.12), in other words, inter- 



ference effects. This can be shown by looking at the den- 
sity of excitation energies, which turns out to be almost 
unchanged by the inclusion of the electron-hole interac- 
tion (see, for example, Fig. 6 in Benedict and Shirley, 
1999). This means that the apparent shift of peak posi- 
tions [about 0.4 eV for the E2 peak in silicon (Albrecht 
et al., 1998a)] is not due to a negative shift of the transi- 
tion energies. 

Of course, this is not true in the case of bound exci- 
tons, which, apart from the above-mentioned insulators, 
includes systems having discrete energy levels, like at- 
oms (Rohlfing and Louie, 2000), molecules and clusters 
(Onida et al., 1995; Rohlfing and Louie, 1998a), and 
other low-dimensional systems. In these cases, excitation 
energies are shifted; moreover, the effect can be very 
strong, since screening is not efficient and the electron 
and hole are localized. 

A good illustration of excitons in low-dimensional sys- 
tems is given by the case of conjugated polymers like 
trans-polyacetylene and poly-phenylene-vinylene 
(PPV). The optical absorption spectra of these mol- 
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ecules have been studied by Rohlfing and Louie (1999a). 
The electron-hole interaction gives rise to exciton bind- 
ing energies of the order of one eV and strongly modi- 
fies the spectral line shapes. Similar results for polymers 
have also been found by van der Horst et al (1999) and 
by Ruini et al. (2002). 

Surfaces can also be considered as low-dimensional 
systems, in particular in certain reconstructions like the 
quasi-one-dimensional chain structures of Si(lll)2xl 
and Ge(lll)(2xl). A strong excitonic effect has been 
found for Si( 111)2x1 in semiempirical (Reining and 
Del Sole, 1991) and ab initio calculations (Rohlfing and 
Louie, 1999b). The excellent agreement found allows 
one to use the theoretical results as reference spectra for 
the interpretation of experiments. For example, an ab 
initio exciton calculation by Rohlfing et al. (2000) for 
Ge(lll)(2xl) has demonstrated how optical differen- 
tial reflectivity spectra can be used to distinguish be- 
tween the two possible isomers of the reconstructed sur- 
face. This distinction was made possible by the fact that 
a quantitative comparison between the calculated and 
experimental spectra is possible when electron-hole ef- 
fects are treated correctly (see Fig. 7). 

Rohlfing and Louie (2000) have recently published an 
extended paper in which their first-principles exciton ab- 
sorption calculations are reviewed. They present the 
theoretical framework, compare different choices of ba- 
sis sets, and discuss the applications to both finite sys- 
tems (namely, He, Ne, and Ar atoms, and SiH 4 and 
larger Si„H„, clusters) and infinite crystals (GaAs, Si, 
LiF, and LiCl). 



The importance of excitonic effects can also be seen in 
electron-energy-loss spectra, linked to the imaginary 
part of the inverse dielectric function. Recently, 
Soininen and Shirley (2000) applied the ab initio exciton 
scheme to the calculation of the dynamic structure fac- 
tor S(q,w) of diamond, LiF, and wurtzite GaN for vari- 
ous momentum transfers q. As in the case of the macro- 
scopic dielectric function, the inclusion of the electron- 
hole interaction in 8 _1 was found to substantially 
redistribute the spectral weight with respect to a GW 
(independent-quasiparticle) calculation, being crucial in 
interpreting the experimental dynamical structure factor 
over a wide range of energies. Essentially, an overall 
shift to lower energies and an enhancement of the inten- 
sities on the low-energy side are found. The latter effect 
can be traced back to the nontrivial excitonic effects in 
the imaginary part of the dielectric function. Concerning 
the overall, essentially rigid shift, one might wish to dis- 
cuss in detail the extent to which the electron-hole at- 
traction and self-energy corrections cancel in the dy- 
namical structure factor. For the plasmon peak of silicon 
at vanishing momentum transfer, it has in fact been 
shown that strong cancellations occur (Olevano and 
Reining, 2001a). In other words, whereas it is clear that, 
at present, absorption spectra of bulk materials most of- 
ten necessitate the inclusion of self-energy and electron- 
hole attraction effects, it is less obvious how much these 
effects improve loss spectra in general and to what ex- 
tent improvements obtained with respect to using Eq. 
(2.8) are essentially due to local field effects. This latter 
point will also be discussed in Sec. VI. D. 

In conclusion, it is important to note that a large va- 
riety of problems has by now been studied via the ab 
initio approach, ranging from absorption to energy-loss 
spectra; from applications to atoms, molecules, and clus- 
ters to applications to bulk materials; and showing exci- 
tonic effects as different as the continuum exciton in 
silicon, the very weakly bound Wannier excitons in 
GaAs, or the strongly bound exciton in LiF. Models exist 
which allow each of these situations to be described, the 
two best known being the hydrogenic Wannier model 
for weakly bound excitons and the Frenkel model for 
strongly bound excitons (see Bassani and Pastori Par- 
ravicini, 1975 for an overview). In fact, certain features 
in the spectra can be very easily predicted by these mod- 
els. In particular, the knowledge of the effective masses 
(band curvatures) at the direct gap, and of the dielectric 
constant, often yield a very good estimate of exciton 
binding energies and allow one to estimate strongly 
bound excitons like those in LiF. However, none of the 
simple models could ever predict an entire spectrum. 
Solid argon, for example, exhibits a hydrogenlike series 
of peaks below the continuum, the first of which is 
Frenkel-like, and the others of which are well described 
by the Wannier model. More refined model approaches 
manage to include both types of bound excitons (Resca 
et al, 1978), but are not meant to predict excitonic ef- 
fects in higher parts of the spectrum, where other details 
of the band structure become important. Instead, the ab 
initio approach has been shown to treat all those fea- 
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tures equally well, and hence to be flexible and predic- 
tive. Of course the agreement with experiment is still 
somewhat limited by the various approximations that 
are commonly used. 

In the following section, several possible choices will 
be discussed. 

4. Different levels of sophistication 

One of the main difficulties one encounters when dis- 
cussing the quality of Bethe-Salpeter results is the fact 
that, as can be seen from Fig. 4, the ab initio approach is 
a three-step method. First, a Kohn-Sham ground-state 
calculation is performed. Second, GW corrections are 
added in order to obtain correct quasiparticle energies. 
Third, the Kohn-Sham and GW results are used in order 
to construct the Bethe-Salpeter equation. At each step, 
approximations are made and one has to be careful in 
order to keep the error bar in the final results small. 

In the following, we shall examine a list of the main 
approximations that are often used. 

(i) Pseudopotentials. All above-cited ab initio calcu- 
lations use pseudopotentials. The use of pseudo 
wave functions for the construction of transition 
matrix elements might, however, introduce some 
error in the results, even when the nonlocality of 
the pseudopotential is correctly taken into 
account. 3-1 The influence of strongly nonlocal 
pseudopotentials on the calculation of transition 
matrix elements has been analyzed by Read and 
Needs (1991), who have found differences of up 
to 10-15 %, depending on the element. See also 
Marini et al. (2001) for an example. 

Arnaud and Alouani (2000, 2001) have calculated 
quasiparticle energies and optical spectra includ- 
ing self-energy and excitonic effects, using the all- 
electron full-potential projector-augmented wave 
method of Blochl (1994). For Si, GaAs, and dia- 
mond their GW band structure agrees with the 
available pseudopotential results to within 1 % for 
GaAs and AlAs, 2% for Si, and 5% for InP (Ar- 
naud and Alouani, 2000). Optical spectra com- 
puted within the Bethe-Salpeter scheme (Arnaud 
and Alouani, 2001) are also in agreement to 
within a few percent with pseudopotential results, 
except for the case of diamond, where they found 
differences in the peak positions of the order of 
0.5 eV. 

(ii) Kohn-Sham-LDA wave functions. In the Bethe- 
Salpeter approach, Pjq P and other matrix ele- 
ments are most often constructed using GW ei- 
genvalues but LDA wave functions. This has been 
justified by the fact that quasiparticle and Kohn- 
Sham orbitals are supposed to be close to each 



j3 Naturally, when using nonlocal pseudopotentials, the cor- 
rect definition of the velocity operator must be adopted, since 
macroscopically wrong results can be obtained otherwise (Van 
Dyke, 1972; Read and Needs, 1991). 



other (Hybertsen and Louie, 1986; Godby et al, 
1986, 1988). However, this is not always the case. 
For example, for finite systems the missing — l/r 
tail of the V xc potential can strongly distort the 
wave functions, and it may be necessary to correct 
their spatial behavior. This has been done by 
Rohlflng and Louie (2000). For the small clusters 
studied, the authors have shown that the wave 
functions of the empty states become clearly more 
delocalized with respect to the LDA Kohn-Sham 
ones. This has a sizable effect on the exciton bind- 
ing energies, which are reduced by as much as 
20%. Note that it would be much more compli- 
cated to include exact quasiparticle wave func- 
tions. In fact, in order to obtain Eq. (4.11), the 
wave functions used for the construction of the 
matrix elements have to derive from a static op- 
erator. Otherwise, they might be nonorthogonal 
for different states (energies), the above basis 
transformation would become more complicated, 
and the calculation could yield different results. 

(iii) Static electron-hole screening. This approximation 
is often justified, since the plasma frequencies of 
the investigated systems are much larger than the 
excitonic binding energies. Moreover, it has been 
found that dynamical effects in the electron-hole 
screening and in the one-particle Green's function 
tend to cancel each other, at least in simple semi- 
conductors (Bechstedt et ah, 1997), which suggests 
that both should be neglected. This might not be 
true in other systems where the exciton binding 
energy is large. Dynamical effects in the Bethe- 
Salpeter equation have been studied for core ex- 
citons by Strinati (1982, 1984), who derived the 
more general effective two-particle equation 34 

2 ^?{E k )Ai c \E x )=E k Al c (E x ). (4.13) 

v'c' 

Note that, as in the case of the one-particle ex- 
cited states, it is the dielectric function that intro- 
duces the energy dependence (and hence finally 
the possibility of complex eigenvalues) into the 
equation: the two-particle excited state may now 
have a finite lifetime, shorter than the electron 
and hole lifetimes. And again, since this comes 
through e _1 , energy is conserved, since one exci- 
tation is decaying by exciting others. Dynamical 
effects have been studied for SiH 4 by Rohlflng 
and Louie (2000); they are found to increase the 
exciton binding energy by about 0.1 eV. 

(iv) RPA electron-hole screening. The kernel W is usu- 
ally evaluated in the RPA. Since the neglect of 
exchange-correlation effects in the screening can 
lead to an underestimate of the dielectric con- 
stant, this approximation might overestimate exci- 



Strinati found that dynamical screening effects significantly 
narrow the core-excitation spectral width, which goes along 
with an increase in the exciton binding energy. 
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ton binding energies. Chang et at (2000) have 
tried to estimate the effect for a-quartz by rescal- 
ing the electron-hole screening self-consistently, 
with the dielectric constant calculated including 
excitonic effects. In fact, they find a shift of the 
main absorption peak towards higher energies, 
closer to the experimental one. A quantitative dis- 
cussion is difficult, since it should include a recal- 
culation of the dielectric matrix for all momentum 
transfers and a consistent inclusion of vertex cor- 
rections (i.e., an additional vertex, because the 
Bethe-Salpeter equation should now be derived 
from 2 = /GWT). 
(v) Neglect of resonant-antiresonant coupling. Gener- 
ally, Eq. (4.12) is used for the calculation of ab- 
sorption spectra. This supposes that transitions at 
positive (resonant part) and negative (antireso- 
nant part) frequencies do not mix (see Appendix 
B.2). The effects of the resonant-antiresonant 
coupling on the exciton binding energy have been 
discussed by Zimmermann (1970) and Del Sole 
and Selloni (1984). Neglecting the coupling can 
bias the results, especially for quantities based on 
the real part of epsilon, like the dielectric constant 
or the electron-energy-loss spectra. This has been 
illustrated by Olevano and Reining (2001a) for 
the case of the plasmon resonance of silicon. 

To summarize, we stress that the Bethe-Salpeter ap- 
proach has up to now yielded results that agree with the 
measured ones within 10% for the peak positions and 
20% for the peak strengths. The error bar is given by the 
sum of several contributions, and it has not yet been 
completely elucidated where the main sources of error 
lie. Well-generated ab initio pseudopotentials should af- 
fect the results by less than —0.1 eV. 35 Also, as men- 
tioned in Sec. IV.A.2, GW transition energies may be off 
by 10-15 %, and one suspects that the quality of the GW 
result dominates the agreement between the experimen- 
tal and the theoretical peak positions. The way the 
screening is taken into account (e.g., with a model di- 
electric function, or in the static RPA) also contributes 
to the error bar by at least another 0.1 eV. As for present 
limitations of the approach, one should mention the ne- 
glect of the dynamical screening of the electron-hole in- 
teraction, which could become relevant when looking at 
the dielectric function at frequencies comparable with 
the plasma frequency, and the complete neglect of the 
lattice vibrations. Finally, there are the limitations due to 
the computational heaviness of the approach, which is in 
part necessarily cumbersome, since it involves four-point 
quantities. This point is, among others, a strong motiva- 
tion to search for alternative approaches, the most 
prominent being the TDDFT approach, to be discussed 
in the next section. 



The choice of the lattice constant — experimental or 
theoretical — can also contribute 50-200 meV. 



V. TIME-DEPENDENT DENSITY-FUNCTIONAL THEORY 

Several reviews on general foundations of TDDFT 
and some applications have appeared recently, including 
those of Gross et al. (1994, 1996), Casida (1995, 1996), 
Dobson, Vignale, and Das (1997), Burke et al. (2002), 
and van Leeuwen (2001). Once more we focus the dis- 
cussion below on the points that are essential for the 
comparisons that are within the scope of this review. In 
the following, we summarize the foundations of TDDFT, 
and its connections with excited-state properties and 
with many-body perturbation theory (a detailed com- 
parison between TDDFT and the Green's-function ap- 
proach will be given in Sec. VI). 

A. Formalism 

The Hohenberg-Kohn-Sham theory as described in 
Sec. II. A is a ground-state theory, and is hence not 
meant for the calculation of electronic excitations. How- 
ever, one can extend the idea of static DFT by analogy 
with classical mechanics: the ground state is determined 
by the minimum of the total energy. When one asks for 
the evolution of the system under the influence of a 
time-dependent external potential, one should search 
the extrema of the quantum-mechanical action, 

A= \ tl dt(V{t)\i^--H{t)\V(t)), (5.1) 
J to dt 

just as in classical mechanics the trajectory of a system is 
determined by the extrema of the action j' t ^dtL(t), 

where L is the Lagrangian. Theorems have now been 
established for time-dependent DFT which are parallel 
to those of static DFT. Many-body effects are included 
in a local exchange-correlation potential, which is now 
time dependent and, as in the case of static DFT, is un- 
known (for more details see Gross et al., 1996). The first 
applications of TDDFT response theory were made be- 
fore the formal development and relied on analogies 
with time-dependent Hartree-Fock theory (Stott and 
Zaremba, 1980; Zangwill and Soven, 1980). In the work 
of Runge and Gross (1984) a theory similar to that of 
Hohenberg, Kohn, and Sham is developed for time- 
dependent potentials in terms of the action functional. 
The first theorem proves a one-to-one mapping between 
time-dependent potentials and time-dependent densities 
(that are v representable); the second proves the 
stationary-action principle. The scheme is very similar to 
the ground-state Kohn-Sham formalism. The proof of 
the first theorem is based directly on the evolution of the 
time-dependent Schrodinger equation from a fixed ini- 
tial many-particle state ^(/ ) = 1 P under the influence 
of a time-dependent potential V exl (t) (required to be 
expandable in a Taylor series around f ). The initial 
state W does not need to be the ground state or any 
other eigenstate of the initial potential. As one does not 
rely on the adiabatic connection in standard zero- 
temperature many-body perturbation theory, the for- 
malism is able to handle external perturbations varying 
rapidly in time. 
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By virtue of the first theorem, the time-dependent 
density determines the external potential uniquely up to 
an additive purely time-dependent function. On the 
other hand, the potential determines the time- 
dependent wave function; therefore the expectation 
value of any quantum-mechanical operator is a unique 
functional of the density. This theorem has been gener- 
alized by van Leeuwen (1999) by showing that a time- 
dependent density obtained from an initial many-body 
system can be, in principle, reproduced by a time- 
dependent external potential in a many-body system 
with different (and possibly null) two-particle interac- 
tion. This mapping between densities and time- 
dependent potentials in TDDFT demonstrates the non- 
interacting v representability of a general many-body 
system. That is, if one takes the two-particle interaction 
of the second system to be zero, this theorem establishes 
that for a given initial state having the proper initial 
density and time derivative, there is a unique potential 
(up to a purely time-dependent function) in a noninter- 
acting system that reproduces the given time-dependent 
density at all times. This result is important for the 
Kohn-Sham formalism of TDDFT (see below). 

Moreover, in addition to their dependence on the 
density, the time-dependent functionals depend on the 
initial state ^ . 36 The time-dependent particle and 
current density can be calculated exactly from the con- 
tinuity equation and the equation of motion of the 

paramagnetic current-density operator j p (r) 
= (l/2i)'2f :=1 [\.S(r-r j ) + S(r-rj)V Ij i, that is, 



dp(r,t) 
dt 

dj(r,t) _ 

dt 



■vjd.o, 



-i(nt)\[L,H(t)]\nt)). 



(5.2) 



(5.3) 



In particular, the current density j(r,r) following from 
the time-dependent Kohn-Sham orbitals and the true in- 
teracting current density may differ only by the curl of 
an arbitrary function (Dobson, Vignale, and Das, 1997). 
The second theorem deals with the variational principle 
of the action functional with the initial condition ^(f ) 
= ^0- From the previous one-to-one mapping between 
time-dependent potentials and densities, the action [Eq. 
(5.1)] is a functional of the density A[p], which must 
have a stationary point at the correct time-dependent 
density. Thus the Euler equation corresponding to the 
extrema of A[p], SA[p]l Sp(r,t) = determines the 
time-dependent density, just as in the Hohenberg-Kohn 
formalism the static ground-state density is given by the 
minimum of the total energy E(SE[p]l <5p(r) = 0). Simi- 
larly, one can define a time-dependent Kohn-Sham 
scheme by introducing a noninteracting system that re- 



j6 In general there are many wave functions leading to the 
same density. Any of them will work properly, because the 
dependence of the effective time-dependent potential is such 
that the interacting density will be reproduced in each case. 



produces the exact interacting density p(r,r). Assuming 
the (demonstrated) v representability of the time- 
dependent densities (van Leeuwen, 1999), one gets the 
following time-dependent Kohn-Sham equations: 



1 , 

2 V +y ' 



elf 



(1,0 



>l'i(r,t) = i—i//i(r,t) 



p(r,t)- 



N 



(5.4) 



(5.5) 



where V elf (r,t) = V H (r,t) + V xc (r,t) + V ext (r,t) is the ef- 
fective time-dependent potential felt by the electrons. It 
consists of the sum of the external time-dependent ap- 
plied field, the time-dependent Hartree term, plus the 
exchange-correlation potential (defined via the equiva- 
lence between the interacting and fictitious noninteract- 
ing systems). The variational principle yields 



V xc (r,t)- 



SAM 
8p{r,t) ' 



(5.6) 



where A Je [p] is the exchange-correlation part of the ac- 
tion functional. 37 

The advantage of the time-dependent Kohn-Sham- 
scheme lies in its computational simplicity compared to 
other quantum-chemical methods such as time- 
dependent Hartree-Fock or configuration interaction 
(Langhoff era/., 1972; Jensen, 1999). Up to here, the 
formalism presented deals with the quantum nature of 
the electrons only. The equations can be generalized to 
treat the quantum-mechanical coupling of the nuclear 
and electronic motions (Kreibich, 2000; Kreibich and 
Gross, 2001), the quantum nature of the electromagnetic 
field, and even superconductors (Gross et al, 1996; 
Kurth era/., 1999). 

B. Excitation energies in TDDFT 

Formally, TDDFT allows the calculation of the 
(bound and unbound) excited-state energies and transi- 
tion probabilities of a many-body system, based on in- 
formation gleaned from an ordinary DFT self-consistent 
calculation. In the time-dependent approach, one stud- 
ies how the system behaves when subject to a time- 
dependent external perturbation. In this case, the sys- 
tem's response is directly related to the iV-particle 
excited states of an A^-particle system, in a manner simi- 
lar to the way the one -particle Green's function is re- 



37 As with the total energy in the static case, here the action 
functional has been decomposed as 

A[p]=2jL 1 /;jdf<^(0l«"(a/»)+sV 2 -y ei( (0l^(0> 

-A H [p]-A xc [p], 

where i/z^t) are the wave functions of the fictitious noninter- 
acting Kohn-Sham system, A H is the time-dependent Hartree 
contribution \ J }didtp(t,t)V H (r,t), and A xc includes all ex- 
change and dynamical correlation effects due to the many- 
body interacting system. A xc is not known and has to be ap- 
proximated. 
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lated to the (N+l)- and (N— 1) -particle excited states 
of the same system (see Sec. III. A). 

Petersilka et al. (1996) have developed a formalism 
for calculating the neutral excitations in finite systems. 
The key formula is a Dyson-like representation of the 
exact linear density response x of an interacting many- 
electron system in terms of the noninteracting Kohn- 
Sham response \o ■ To prove this relation, one can start 
from the physical definition of the retarded linear- 
response function, 



SV ext (r',t') 



(5.7) 



which measures the degree to which the density re- 
sponds to first order in the external potential. Equiva- 
lently, the linear response ^ of the fictitious Kohn-Sham 
system that can be described in terms of the Kohn-Sham 
orbitals is given by 



Sp(r,t) 



SV ett (r',t') 



(5.8) 



where V eit =V ext +V H +V xc . Using Sp/SV ext = (Sp/SV ef[ ) 
x (^eff/'5^ / £'.Tr)=A'o^eff/^e.rt one obtains 

SV eii (r,t) 



8V ext (t'f) 



S(r-r')S(t-t') 



S(t-t") 



+f xc (r,t,t",t") 



Xx(*",t",r',t')dr"dt", (5.9) 
where the time-dependent exchange-correlation kernel 
SV xc [p(r,t)] 



f xc (t,t,r',t')-- 



Sp(r',t') 



(5.10) 



has been introduced. Now it is straightforward to get the 
Dyson-like equation 

x(r,r',w) = xo(r,r',w) 



+ J dr 1 dr 2 Xo(r,r 1 ,(o)K(r 1 ,r 2 ,(o)x(r 2 ,r',(o), 

(5.11) 

which has to be solved iteratively and where the kernel 
K has been introduced as 



1 



r 2 



+/*c0l, r 2>«)- 



(5.12) 



Equation (5.11) corresponds to Eq. (2.17) of the In- 
troduction. In fact, Eq. (5.11) and Eq. (5.10) can be ob- 
tained directly, as discussed in the Introduction, for the 
case in which the total potential in the time-dependent 
Hamiltonian is V H + V xc . As also pointed out there, due 
to the density-only dependence of these potentials the S 
functions that appear in the derivation are such that only 
two-point functions are involved from the beginning. 
Equation (5.11) is hence the contraction of Eq. (2.20). 
Of course, the nontrivial point is that the response func- 



tion x calculated from the fictitious noninteracting sys- 
tem using Eq. (5.11) is equal to the true response func- 
tion of the interacting system. 

This scheme provides an exact representation of the 
full interacting linear density response as 



p 1 (r,co)= dr'xo(r,r';u)V eif (r',co) 



- J dr' x (r,r';co)V exl (r',w). 

In other words, the exact linear density response of an 
interacting system can be written as the linear density 
response of a noninteracting system to the effective per- 
turbation V eS . The exact time-dependent exchange- 
correlation kernel is of course unknown, and practical 
calculations must rely on some approximation. The most 
commonly used one, due to its simplicity, is the adiabatic 
local-density approximation, also called the time- 
dependent LDA (TDLDA), in which / re (r 1 ,r 2 ,w) is ap- 
proximated by the (co-independent) functional deriva- 
tive of the LDA exchange-correlation potential: 



n 



TDLDA 



{r l ,r 2 )=c?(r 1 -r 2 ) 



c vy^(p(r 1 ),r 1 ) 
dp{r x ) 



(5.13) 



Apart from this approximation for f xc , another approxi- 
mation has to be made in practical calculations: the 
static Kohn-Sham orbitals and eigenvalues used to con- 
struct xo are in fact calculated with an approximate 
exchange-correlation potential V xc , typically the same 
as the one used in ground-state calculations. 

Formally inverting Eq. (5.11), one obtains a compact 
two-point matrix equation, 

X (co) = [l- xoMKico^xoico), (5.14) 

where the problem of finding the excited-state energies 
of an interacting system (poles of x) has been mapped to 
searching the values of w for which the operator R(to) 
= 1— xo(w)K(a)) is not invertible. In fact, x has poles at 
the true excitation energies a> = O, • , while xo has poles at 
the Kohn-Sham eigenvalue differences. Hence the sin- 
gularities of x must be canceled by the zeros of R(to) 
[and those of xo by the zeros of R~ 1 (co)]. The true ex- 
citation energies ilj can be characterized as those fre- 
quencies where the eigenvalues of R vanish. In other 
words, the energies of the resulting electron-hole excita- 
tions will be renormalized with respect to those of the 
noninteracting electron-hole pairs. This formalism pro- 
vides a convenient starting point for calculating the ex- 
citation spectrum. 

Of course, for the practical solution of Eq. (5.14) all 
points discussed in the Introduction apply. In particular, 
when adding a finite small imaginary part to the fre- 
quency, the two-point equation can be inverted for each 
frequency in order to obtain the spectrum. Alternatively, 
Eq. (5.14) can be transformed into an effective eigen- 
value problem, which implies that one is now working 
with four-point quantities. This second procedure is 
equivalent to solving the generalized four-point eigen- 
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value problem 4 i?(w)|\) = (see Appendix, B.2, in par- 
ticular 4 R = Tl^ 1 = [l— 4 xoK]) and yields the equation 

„ K ( ni n 2 ),(n 3 n 4 )(u) _ 

(5.15) 

where the kernel ^(n 1 « 2 ),(n 3 « 4 ) 

= (®n in2 \K(aj)\<i> n3n4 ), with $ ni „ 2 (r):=^ Bi (r)^* 2 (r). 

The index n embodies the spin as well as orbital degrees 
of freedom. 38 

Looking at the Coulomb contribution to K, one rec- 
ognizes the effective exchange interaction between elec- 
tron and hole (electron-hole exchange), which also ap- 
pears in the Bethe-Salpeter equation (4.11) (for \ an d v 
instead of P and v). The diagonal element reads 

J J d ri d r2 ^* in2 ( ri )^—^^ nin2 ( r2 ), (5.16) 

whereas the f xc contribution within the simple TDLDA 
introduces a local and static attractive electron-hole in- 
teraction: 

f f ^V X c( r l) 

J J dndrafc^ri^ii-ra) — — — <I»„ 1 „ 2 (r 2 ) 

= J drp Bl (r)^^ P „ 2 (r). (5.17) 

One can see the effects of these two different contribu- 
tions on the calculated optical spectrum of a finite sys- 
tem, e.g., a small cluster: the Coulomb part shifts the 
independent electron spectrum to high energies, 
whereas the exchange-correlation brings it partially back 
(this can to a certain extent be compared with the exci- 
tonic corrections computed in the framework of the 
many-body theory, as will be discussed in Sec. VI). 
These effects are clearly seen in Fig. 8, where the optical 
spectrum calculated within TDLDA for three different 
clusters is reported. 

Equation (5.15) can be simplified by an expansion 
around one particular Kohn-Sham transition 1— >2, i.e., 
by calculating the transition energy a> = Cl in first-order 
perturbation theory (single-pole approximation). Fol- 
lowing Petersilka et al. (1996) and considering explicitly 
the spin degrees of freedom, 39 one obtains from Eq. 
(5.15) 

fl — W 12 +Re[i^ 1T 2t,lt2T( w 12) :!: -^lT2|,lT2|( w 12)]- 

(5.18) 



The transition-space representation is based on similarities 
with the quantum-chemistry time-dependent Hartree-Fock ap- 
proach (Langhoff et al, 1972; Casida, 1995). However, a simple 
two-point matrix problem has been transformed into a more 
complex four-point representation. Note also that Eq. (5.15) 
can be derived from the condition det (H 2p — E x )=0 with H 2p 
of the form of Eq. (B20), and K of the form defined in Eq. 
(2.15) without the term involving w. 

39 In the present section we treat spin effects explicitly, as is 
natural when considering atomic transitions. 




Energy (eV) 

FIG. 8. Calculated optical absorption spectrum (given as the 
dipolar strength function) for several clusters (Marques et al. , 
2001), compared to experimental data (Itoh et al., 1986; Wang 
et al, 1990): solid black line, full TDLDA calculations; dashed 
line, independent Kohn-Sham particle approximation (xo)l 
dotted line, an RPA calculation, i.e., putting f xc = in Eq. 
(5.12); gray line, the experimental results. In the inset of Na g 
are the results for the sodium dimer: heavy solid line, full 
TDLDA calculation; gray line, experimental results (Sinha 
et al, 1949), but now the dashed line (almost indistinguishable 
from previous one) is for a calculation using an exact-exchange 
functional (Marques et al, 2001). As discussed in the text the 
kernel includes an effective attractive part (electron-hole at- 
traction), which reduces the Coulomb term (electron-hole ex- 
change). 



Note that the noninteracting Kohn-Sham response func- 
tion is diagonal in the spin variables and exhibits poles 
at the Kohn-Sham energy differences corresponding to 
noninteracting electron-hole excitations within the same 
spin space. The mixing of spin channels comes into play 
simply by the spin-dependent exchange-correlation ker- 
nel, and the magnetization response naturally involves 
spin-flip processes. In order to make more explicit the 
fact that the approximation in Eq. (5.18) embodies the 
spin-multiplet structure of the excitation spectrum of 
otherwise spin-unpolarized ground states, one can re- 
write the f xc kernel in terms of the two independent 
combinations of the spin components of the kernel: f xc 
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=0.5(/« n +/« u ) and/2 t = 0.5(/ xcTT -A. cT1 ). 40 Thus the 
singlet and triplet solutions of Eq. (5.18) are (note again 
the similarity to the Bethe-Salpeter equation concerning 
the contribution of v) 

o>sin g iet= w 12 +2Re j t3?r 1 dr 2 3>f 2 (r 1 ) 

X ( \ t - T \ + Ac( r l ' r 2 ^12) ) * 12(^2) 

and 

(0 triplet = w 12 

,wi 2 )'I ) i 2 (r 2 ). 

Promising results are obtained for the lowest excitation 
energies of atoms and molecules by using different ap- 
proximations for V xc and f xc (Petersilka etal., 1996; 
March etal, 1999; Casida etal., 2000; Petersilka etal, 
2000). 

Results for finite systems beyond the single-pole ap- 
proximation which have been particularly discussed are 
those of TDLDA (Rubio etal, 1996; Yabana and 
Bertsch, 1996, 1999a, 1999b; Vasiliev etal, 1999, 2002; 
Casida et ah, 2000; Marques et al, 2001) and the results 
obtained within the OEP scheme 41 (Gorling, 1999; 
Gross etal, 1996; Petersilka etal, 1996, 2000; Grabo 
etal, 2000a, 2000b). The latter approach gives an 
exchange-correlation potential with the correct —llr be- 
havior at long distances (Casida and Salahub, 2000). The 
optimized effective potential yields a uniform shift of 
the energies of transitions to Rydberg states (Petersilka 
etal, 2000; Stener etal, 2001) that mimics to some ex- 
tent the results obtained using the exact V xc available 
for some atoms (Petersilka et al, 2000). From this work 
the importance of a good description of the static 
exchange-correlation potential is clear, as is also the fact 
that only ground-state quantities are needed in the cal- 
culation of excitations. V xc is especially important for 
the higher atomic excited states, which are almost uni- 
formly shifted from the true excitation energies (the 
shift can be related to the difference between the ioniza- 
tion potential and the highest occupied orbital, which in 
exact DFT should be zero). Further developments of 
time-dependent functionals to treat problems involving 
excited-state dissociation (Cai et al. , 2000; Aryasetiawan 
et al. , 2002a) and autoionization resonances (Stener 
etal, 2001) are needed. In the latter case, even if the 
atomic photoionization cross sections are well described 
using the exact V xc and the TDLDA kernel, this is not 



40 The f xc part of the kernel describes exchange-correlation 
processes in the Kohn-Sham system related to the linear re- 
sponse of the frequency-dependent magnetization density, 
whereas f xc is related to the frequency-dependent density. 

41 The calculations in the OEP scheme handle exchange ex- 
actly, and correlations are treated in a local or gradient- 
corrected adiabatic functional. 



the case for the autoionization resonances that turn out 
to be very sensitive to the particular choice of f xc . 

In extended systems with delocalized states one can 
directly solve Eq. (5.14) in the complex plane. [The 
secular Eq. (5.15) becomes a four-point integral equa- 
tion and is computationally not advantageous]. The ze- 
ros of R provide both the excitations as well as the cor- 
responding lifetimes. This technique has been 
successfully used to describe the plasmon dispersion of a 
homogeneous electron gas (Tatarczyk et al, 2001). How- 
ever, the TDLDA has only limited success in describing 
excitations in extended systems, and most likely both 
nonlocality in space and time are needed in order to get 
this renormalization. The successes and failures of 
TDLDA in finite and infinite systems, respectively, are 
not yet fully understood. However, it appears that im- 
provements might more easily be found through an im- 
proved V xc in the case of finite, and an improved f xc in 
the case of infinite, systems. In any case, f xc is a quantity 
of interest which (i) has been discussed in less detail 
than V xc and (ii) is necessarily more complicated, since 
it is its functional derivative. Therefore we look more 
closely at f xc in the following section. 

C. The exchange-correlation kernel f xc 

In the previous section we introduced the exchange- 
correlation kernel f xc to take into account all the dy- 
namical exchange and correlation effects in the response 
of a system to an external perturbing potential. An exact 
representation of f xc in terms of the response functions 
is obtained directly from Eq. (5.11): 

f xc (r,t;r',t') = x 1 (r,t;r',t')-x~ l (r,t;r',t') 



|r— r I 

Note that for finite systems the frequency-dependent re- 
sponse operators can be noninvertible at isolated fre- 
quencies (isolated real poles). However, this is no longer 
the case for infinite bulk systems. In Appendix A we 
provide a summary of the known exact properties of f xc 
that can be used as constraints to build new approxima- 
tions to the exchange-correlation kernel. 

As a consequence of causality, response functions 
must be zero for t'>t. Therefore the f xc kernel cannot 
be symmetric under the exchange of (r,t) and (r' ,t'). 
Since, on the other hand, f xc is the functional derivative 
of the exchange -correlation potential, one concludes 
that either the exact V xc (r,t) cannot be a functional de- 
rivative of the A xc functional (i.e., f xc is not a second 
functional derivative; see Gross etal, 1996; van Leeu- 
wen, 1998, 1999), or there is a contradiction with the 
stationary-action principle that we have described as a 
basic ingredient of time-dependent density-functional 
theory. This problem applies to all twice-differentiable 
action functionals defined with respect to the physical 
time. It appears when applying the variational principle 
to the action in Eq. (5.1). In fact, from the Runge-Gross 
(1984) theorem, the wave function is determined up to a 
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time-dependent phase factor; this makes the action not 
uniquely defined unless one fixes this phase factor (in 
particular it can be chosen such that A[ 1 I'] = 0). There- 
fore the action is not a density functional as the density 
by itself does not fix the phase of the wave function (van 
Leeuwen, 2001). This apparent contradiction has been 
formally resolved by van Leeuwen (1998, 2001) by de- 
fining a new action functional 42 that properly incorpo- 
rates causality effects, and it is not made stationary but 
rather used as a generating function for the density and 
response functions as in statistical mechanics. 

One of the most widely used approximations to time- 
dependent phenomena is the adiabatic LDA or time- 
dependent LDA (TDLDA) in which the static LDA 
functional is used for the dynamical properties. This 
means that the f xc kernel is a contact function in time 
and space [Eq. (5.13)]. Thus/ Vt . is not frequency depen- 
dent at all. TDLDA gives rather accurate results for sys- 
tems with rapidly varying densities such as atoms, sur- 
faces, and clusters. 43 Gross and Kohn (1985) have 
presented an extension of the LDA approximation to 
include dynamical effects in the f xc kernel, and we refer 
the reader to their work for the details of 
parametrization 44 The idea is to use the homogeneous 
electron-gas kernel /J° m (|r— r'|,&>) and make the as- 
sumption that the linear-induced density is a slowly 
varying function (as in the traditional static LDA). This 
amounts to the approximation / vc .(r,r',w) 
=f h x ° c m ( P (r),w), where f h x ° m is the q = Fourier com- 
ponent ofj* x ° m (p(r),\r'-r"\). 

Other kernels proposed in the literature are the fol- 
lowing (spin variables are omitted for simplicity): 

• Petersilka, Grossman, and Gross (PGG) kernel. This 
was derived by Petersilka et al. (1996) in the context 
of exact exchange time-dependent optimized effective 
potentials and kernels (Gross et al, 1996; Gorling, 
1998) (x-TDOEP), and it is equivalent to the so-called 
Slater approximation in Hartree-Fock calculations 
(Langhoff et al, 1972). It is a frequency-independent 
kernel that in real space reads 



42 The new action is defined using the time contour method of 
Keldysh (1965) in which the physical time is parametrized in 
terms of a parameter called pseudotime. By construction, the 
response functions obtained as higher-order derivatives of the 
action functional are symmetrical in the Keldysh pseudotime. 
A back-transformation to the physical time directly provides 
the desired causal response functions (see van Leeuwen, 2001 
for mathematical details of this technique). 

43 See, for example, (Stott and Zaremba, 1980; Zangwill and 
Soven, 1980; Dobson et al, 1988; Tsuei et al, 1990; Casida, 
1995; Petersilka et al, 1996; Rubio et al, 1996, 1997; Liebsch, 
1997; Vasiliev et al, 1999, 2002). 

44 Note that the lower the density the larger the frequency 
dependence of the kernel. The parametrization can be ex- 
tended to nonvanishing q and to include spin polarization. 
However, it fails to reproduce some exact relations, such as the 
harmonic-potential theorem (Dobson, 1994). 
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2 |r-r'| p(r)p(r') 
Note that the general expression of the TDOEP ker- 
nel (similar to time-dependent Hartree-Fock) is fre- 
quency dependent even at the exchange-only level 
(Petersilka et al., 1996, 1998), a feature that is not ac- 
counted for in the PGG approximation. However for 
a two-electron system at the exchange-only level the 
TDOEP and PGG kernels are equivalent, explicitly 
showing that the dynamical part of the kernel stems 
from processes involving multiple valence transitions. 
TDOEP-SIC kernel. This is an attempt to improve 
over both TDLDA and exact exchange (TDOEP) by 
correcting the TDLDA for the self-interaction error 
( SIC = self-interaction-corrected). This correction 
does not affect antiparallel spin contributions. In this 
case, the kernel reads 



f" 



TDOEP-SIC 



(r,r')=/™^(r,r') 



, 2 f k \M*)*l>t(*')\ 2 

1 k 

2 P (r)p(r') 



r— r 



where p k is the density of the orbital k. This expres- 
sion is exact in the one-electron case and, for more 
electrons, it corrects the spurious self -inter action of 
parallel spin, but keeps the antiparallel contributions 
as in TDLDA. This functional is nevertheless ill de- 
fined, since it is not invariant upon a unitary transfor- 
mation of the Kohn-Sham wave functions. 
BPG kernel. Burke et al. (2002) proposed a hybrid 
method to improve the excitation spectra of small at- 
oms by combining the previous PGG expression for 
symmetric spin orientations and the TDLDA for an- 
tisymmetric spin orientations. For the case of a homo- 
geneous electron gas it reads fz (q) 
= 0.5*^(q)+/™ T ^(q)]. 

CDOP kernel. Corradini et al. (1998) gave a param- 
etrization of the quantum Monte Carlo data of Mo- 
roni et al. (1995) for the homogeneous electron gas 
that satisfies the theoretically known limits for small 
and large q. This kernel has an analytical space Fou- 
rier transform that simplifies its implementation in 
standard first-principles techniques. 
RA kernel. Richardson and Ashcroft (1994) proposed 
a dynamical parametrization of the kernel based on a 
summation of self-energy and fluctuation terms in the 
diagrammatic expansion of the polarization function 
for the homogeneous electron gas. It is constructed to 
satisfy many known exact conditions (static and dy- 
namic). This parametrization is assumed to be very 
close to the exact dynamical exchange-correlation ker- 
nel of the homogeneous electron gas. A corrected ver- 
sion of the RA parametrization can be found in Lein 
et al. (2000). 
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TP kernel. Tokatly and Pankratov (2001), based on a 
many-body diagrammatic expansion, have derived an 
expression for the f xc kernel in terms of the regular 
part Xr °f the Kohn-Sham response function xo with 
respect to a given resonance frequency w ;; : 
f xc {r,r' ,a>ij) 



\ 



dx x dx 2 x r l (r,r 1 )<S> ij (r 1 )<S>f j (r z )x r Hr 2 y) 



where A (7 =(<t> /; |/ YC (w !7 )|$ ;; ) [<J> ;J - was defined after 
Eq. (5.15)]. It is shown that the spatial nonlocality is 
strongly frequency dependent and diverges for the 
case of infinite systems at the excitation energies. This 
result is similar in form to the correction obtained by 
Gonze and Scheffler (1999) and it is also equivalent to 
the Gorling-Levy perturbation theory (Gorling and 
Levy, 1994). 

• RORO kernel. This is a static kernel derived by Rein- 
ing et al. (2002) from a direct comparison between the 
time-dependent density-functional and Bethe- 
Salpeter equations. In Appendix C we develop this 
kernel in greater detail by formally comparing the 
Bethe-Salpeter equation and TDDFT. The kernel 
consists of a contribution stemming from the energy 
shift between Kohn-Sham and GW eigenvalues, and a 
second one describing the electron-hole interaction. 
One can absorb the first, positive contribution into an 
energy shift of the starting Important excitonic 
effects are then obtained by using only the static long- 
range term A/ JCC (q,G,G') = - 5 G G -a/|q+G| 2 . 

Some of these kernels have been tested in two limiting 
cases, for helium and beryllium atoms (Petersilka et al. , 
2000) and for the correlation energy (Lein et al, 2000) 
and plasmon dispersion (Tatarczyk et al, 2001) of a ho- 
mogeneous electron gas. In the case of atoms the de- 
tailed form of the V xc potential is the crucial part in 
getting the absolute position of most excitation energies, 
and the TDLDA kernel provides reasonable results. 
The results are marginally improved using more compli- 
cated kernels, always keeping the "exact" V xc fixed. In 
fact the results obtained by Petersilka et al. (1998, 2000) 
for the helium and beryllium atom using different ker- 
nels indicate that the influence of the f xc on the calcu- 
lated spectra is less important than the choice of a good 
V xc potential. However, this is no longer true for the 
lower excitation energies of beryllium and for singlet- 
triplet splittings, where the effects coming from V xc can- 
cel. Furthermore, Aryasetiawan et al. (2002a) looked at 
the singlet excitation of the H 2 molecule as a function of 
the internuclear distance. The results are summarized in 
Fig. 9, where it can be seen that the simple TDLDA is 
good for intermediate distances only. The inclusion of 
spin dependence improves the large-distance results. By 
comparing with exact results of a simple two-site Hub- 
bard model it was found that, indeed, f xc should be 
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FIG. 9. The excitation energy A.E in the hydrogen molecule 
for the transition from the state to the state as a 
function of the nuclear separation d. LDA stands for the ei- 
genvalue difference in LDA, whereas TDLDA and TDSLDA 
correspond to the calculation of the excitation energy within 
time-dependent LDA without and with spin polarization, re- 
spectively. The full line corresponds to the exact results, pro- 
vided for comparison (Aryasetiawan et al, 2002). 

strongly nonlocal in space and has an important energy 
dependence. The nonlocality has been clearly shown by 
Baerends (2001). He proposed an orbital-dependent 
exchange-correlation potential for the H 2 molecule that 
yields the exact dissociation regime within the Kohn- 
Sham formalism. 

Similar conclusions were obtained by Marques et al. 
(2001) in a study of the optical absorption spectrum of 
small sodium and silicon hydrogenated clusters. 45 More- 
over, it is observed that exchange-correlation kernels fit- 
ted to reproduce atomic properties perform poorly in 
the case of an extended electron gas, due mainly to in- 
correct behavior at long wavelengths. This limit is well 
reproduced in the TDLDA but not in the PGG kernel. 
Furthermore, the static CDOP kernel gives results of 
similar quality to those of the more elaborated dynami- 
cal RA kernel. This indicates that the frequency depen- 
dence of the kernel is of little importance in providing 
good total correlation energies and plasmon dispersion. 
Note that TDLDA does not perform too badly in this 
last case (Larson et al, 1996). 

Finally, Reining et al. (2002) have tested the above- 
mentioned long-range contribution to the static RORO 
kernel by performing a TDDFT calculation for bulk sili- 
con in the following way. First, they determined the 
LDA electronic structure. Second, they constructed xo , 
but with the eigenvalues shifted to the GW ones, in or- 
der to simulate the first part of the kernel as explained 
above. Third, they used A/ xc (r,r') = — a/(4ir\r— r'|), 
with the empirical value a = 0.2. The result of the 
TDDFT-RORO calculation is shown in Fig. 10. The dots 
are the experimental results for the absorption spectrum 
measured by Lautenschlager et al. (1987). The dotted 
and dot-dashed curves are the results of a RPA and a 



45 All the results for the optical spectrum of small Na clusters 
are very similar, regardless of the exchange-correlation poten- 
tial used. In contrast, hydrogenated silicon clusters show that a 
much better agreement with diffusion quantum Monte Carlo 
calculations (Grossman et al., 2001; Porter et al., 2001) and 
with experiments can be obtained when the exact exchange 
potential is used. 
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co[eV] 

FIG. 10. Optical spectrum [Im(e M )] for silicon: 9, experiment 
(Lautenschlager et al., 1987); dotted curve, RPA; dot-dashed 
curve, TDLDA; long-dashed curve, Bethe-Salpeter equation; 
continuous curve, TDDFT using RORO kernel (Reining et al, 
2002; see text for details). 

TDLDA calculation. The well-known discrepancies with 
experiment are found. The solid curve is the result of the 
approximate TDDFT-RORO calculation, visibly in close 
agreement with the Bethe-Salpeter result (dashed) and 
with experiment. 

Similarly good agreement has been found for the op- 
tical spectra of other semiconductors (Botti et al, 2002). 
It turns out that this static long-range contribution to the 
kernel is sufficient to reproduce strong continuum exci- 
tonic effects in semiconductors. Promising results were 
also obtained for the optical spectrum of some simple 
semiconductors by Kim et al. (2002) using the exact ex- 
change Kohn-Sham formalism together with a TDLDA 
f xc , and by de Boeij et al. (2001) using a polarization- 
dependent functional within the current-density func- 
tional proposed by Vignale and Kohn (1996). This pro- 
cedure shows the influence of macroscopic electric fields 
in the response function. However, none of these ap- 
proaches has up to now managed to describe bound ex- 
citons, and one should certainly at least go beyond the 
simple — 1/q 2 approximation for the RORO kernel in 
order to describe such effects. 

One can conclude that the use of TDDFT for the cal- 
culation of neutral excitations is promising: in the low- 
energy range of the absorption spectra of clusters, 
TDDFT, even in the adiabatic local-density approxima- 
tion, significantly corrects peak positions with respect to 
those found when Kohn-Sham eigenvalue differences 
are interpreted as excitation energies. Typically, the er- 
ror in excitation energy reduces from an amount of the 
order of 50% to an order of magnitude of 10% (often, 
the main correction comes from v or v, i.e., the variation 
of the Hartree potential). Also, valence plasmons in sol- 



ids are well described in TDLDA for small, momentum 
exchange, yielding plasmon frequencies that agree with 
the experimental ones within 0.5-2 eV, and good line 
shapes. However, the optical absorption spectra of solids 
are not improved by TDLDA with respect to the sum of 
Kohn-Sham transitions (see also Sec. VI.D.l). More pre- 
cisely, calculations performed on many different semi- 
conductors (Gavrilenko and Bechstedt, 1996; Kootstra 
et al. , 2000) demonstrate that TDLDA spectra are, on 
average, redshifted by about 40% of the gap. Even when 
the gap is corrected with a scissor operator [as is com- 
monly done for LDA calculations (Levine and Allan, 
1989)], significant discrepancies with experiment remain: 
typically, in zinc-blende semiconductors the strength of 
the E x peak is underestimated by 0-40 %, that of E 2 is 
overestimated by 30-100 %, and the absorption onset is 
overestimated by 0.5 eV (e.g., in GaSb, InAs, ZnTe) or 
even by 0.8 eV (in CdTe and InSb). It is clear that better 
exchange-correlation potentials V xc and kernels f xc must 
be found in order to make TDDFT become the ap- 
proach for the calculation of absorption spectra, just as 
static DFT is today the predominant method for ab ini- 
tio calculations of ground-state properties. In fact, f xc 
must be a strongly nonlocal functional of the density, 
and, in principle, frequency dependent. It seems to be 
crucial to include in f xc at least a contribution going as 
— 1/q 2 for >0, in order to reproduce continuum exci- 
tonic effects in the absorption spectra of infinite systems. 
See also Appendix C. 

VI. TDDFT VERSUS BETHE-SALPETER 

TDDFT and the Green 's-f unction approach represent 
two complementary ways to calculate a dynamical di- 
electric matrix and the spectra derived from it. The ulti- 
mate goal should of course be either to decide whether 
one of the two approaches is clearly superior in a given 
situation, or whether it is possible to combine the advan- 
tages of both in order to design a reliable and efficient 
way to calculate electronic spectra for a broad range of 
applications. It is therefore useful to compare the two 
approaches on different levels: from the mathematical 
point of view, by trying to work out which effects the 
different approximations might cause, and by comparing 
results for different systems and types of spectroscopies. 
It should not be forgotten that in principle both ap- 
proaches are exact, and failures can only be explained 
by the unavoidable approximations, or by an application 
that is in principle not adequate for a particular problem 
(as in the case of static DFT and the interpretation of 
Kohn-Sham eigenvalues as electron addition or removal 
energies). 

A. The equations 

Let us first examine the structure of the equations. It 
should be kept in mind that, whether absorption or loss 
spectra are to be calculated, and whether TDDFT or the 
Bethe-Salpeter approach is used, the equation to be 
solved is of the form S = P]q P + Pjq P KS [see Eqs. (4.10) 
and (5.12)]. Here S may be either x or P (describing loss 
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and absorption spectra, respectively), depending on 
whether the kernel K contains the full Coulomb interac- 
tion v or the truncated v. Moreover, K has a contribu- 
tion F containing the exchange-correlation effects. 
There are two important differences between the 
TDDFT and the Bethe-Salpeter approaches. First, -P/gp 
is either the Kohn-Sham (^ ) or the quasiparticle 
independent-particle response. The second difference 
concerns F. In particular, since the Bethe-Salpeter ap- 
proach derives from Green's functions (i.e., density- 
matrix-like quantities), it naturally leads to a four-point 
function F, whereas the TDDFT equation, which is 
dealing with density-only potentials, can be contracted 
and yields a two-point equation, since the delta func- 
tions in the v or v part of the kernel [Eq. (4.10)] are such 
that one can contract Eq. (4.9) and get the two-point S 
without solving the four-point equation first [this corre- 
sponds to using the f xc defined in the previous section, 
Eq. (5.10)]. 

The Bethe-Salpeter equation with the full kernel [Eq. 
(4.10)] cannot be contracted, because of the fact that the 
S functions of the exchange-correlation part are connect- 
ing different indices from those of the local field part. If 
one supposes for a moment that the self-energy could be 
approximated by the DFT exchange-correlation poten- 
tial, for the calculation of one-quasiparticle eigenvalues 
but also for the functional derivative leading to the ker- 
nel, one would obtain 



<SG(r 3 ,r 4 ) 



5(rj-r 2 ) dr 5 



SV xc ( ri ) Sp(r 5 ) 



Sp(r 5 ) <5G(r 3 ,r 4 ) 



z fxc(*i ,r 3 ) <?(!-! -r 2 )<5(r 3 -r 4 ), 



(6.1) 



the same structure of connecting S functions as in the 
Hartree (i.e., v) part. Again, one can contract the equa- 
tion and obtain the TDDFT two-point form. This is of 
course true for any form of the exchange-correlation po- 
tential that is purely density dependent and local, 
whereas for the nonlocal and Green's-function- 
dependent potentials of many-body perturbation theory 
the equations are necessarily four-point ones. Of course 
this reasoning should not cause the misunderstanding 
that TDDFT is an approximation to the Bethe-Salpeter 
equation; besides the fact that it is a true and in principle 
exact alternative, there are also subtleties linked to the 
use of time-ordered and retarded quantities which pro- 
hibit an easy switching between the two. However, the 
discussion of the structure of the equations still remains 
valid. 

Equations (4.9) and (5.11) are not the only way to 
write the Bethe-Salpeter and TDDFT schemes. In fact, 
in the previous sections we have already seen that it is 
often convenient to transform these equations to an ef- 
fective eigenvalue equation, namely, Eq. (4.11) for the 
Bethe-Salpeter equation, and Eq. (5.15), multiplied by 
(co—o)j k ), for TDDFT. It should be pointed out that in 
this case both equations have become four-point equa- 



tions, because they have been obtained by a basis trans- 
formation to a basis of pairs of occupied and empty 
states. The two equations now look exactly identical, the 
only differences being the meaning of the one- 
(quasi)particle eigenvalues and the content of the ker- 
nel. Clearly, this way of rewriting the equations is con- 
venient if (a) the single diagonalization to be performed 
is less onerous than the inversion of x Ior many frequen- 
cies (which necessarily assumes that the kernel does not 
depend on frequency), and/or (b) if the basis of pairs of 
states which must be considered is smaller than a real- 
space (or reciprocal-space) basis. This is naturally true 
for the Bethe-Salpeter equation, where those basis sets 
would also be quadratic. In the case of TDDFT, the qua- 
dratic basis of states must be compared to the linear 
basis in real or reciprocal space. In that case, the choice 
of the four-point form can only be convenient for small 
finite systems with well-spaced energy levels and in 
which only a limited part of the spectrum is needed. A 
larger energy range would require a higher number of 
states, increasing considerably the computational effort 
in this four-point approach. This explains why one finds 
this form mostly in applications of quantum chemists. 

Writing the equation in the basis of transitions has the 
additional drawback that the full Hamiltonian [Eq. 
(B20)], i.e., that including resonant, antiresonant, and 
coupling contributions, is non-Hermitian. It has been 
shown that this problem can be avoided by transforming 
the equation into a simple quadratic one of the size of 
the resonant contribution only (Casida, 1995; Bauern- 
schmitt and Ahlrichs, 1996): 

(R-C)(R + C)a = co 2 a, (6.2) 

where R and C are the resonant and coupling matrices, 
respectively, and a is a linear combination of the two 
corresponding parts of the full eigenvector A k (see Ap- 
pendix B.2 for details). In the case of TDDFT and for 
real wave functions, R vcv , c , = (o uc S vv ,S cc , + K vcv , c , and 
C vcv i c > = K vcv , c , . Then R — C is diagonal, and the qua- 
dratic equation can be simplified; it reads (in its symme- 
trized form) 



[coj k SjiS km + 2\lfj k (o jk K jk Jm (co) \lfi m coi m ]c 



■■ co 2 c jk . 



(6.3) 

where all the spin degrees are embodied in the i,j in- 
dexes, and fjk = fj~fk ■ Here one is explicitly taking into 
account transitions of the type I— >m where / is occupied 
and m is unoccupied. In the Bethe-Salpeter equation, 
the electron-hole attraction term does not show the sym- 
metry C VCV , C , = K VCV , C , , and Eq. (6.2) cannot be simpli- 
fied. This can be a considerable problem in applications 
like the calculation of electron-energy-loss spectra, 
where the resonant-antiresonant coupling cannot be ne- 
glected, since in that case the non-Hermitian Hamil- 
tonian [or, equivalently, Eq. (6.2)], do not allow the ap- 
plication of fast inversion techniques like the Haydock 
recursion method (Haydock, 1980). 

The third main representation of the TDDFT equa- 
tion is the explicitly time-dependent one. In fact, the 
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spectra we are interested in are due to the response of 
the system to a time-dependent external field. It is hence 
possible to extract the desired information from a calcu- 
lation of the time evolution of the density when such a 
perturbation is applied to the system. This corresponds 
to going one step back in the derivation of the equa- 
tions, as was done, for example, in the introduction to 
this review. Since in TDDFT the density is easily con- 
structed from the time-dependent wave functions, it is 
sufficient to solve the time-dependent Schrodinger equa- 
tion for ifj and construct the density at every step. This is 
particularly easy when approximations for the potential 
are chosen that do not contain memory effects, which 
means that not only are the requirements of computer 
time modest but the storage requirements are low. This 
method thus uses a direct solution of the time- 
dependent single -electron Schrodinger equation for the 
occupied states, 



i-^^ = H KS (t)^(r,t) (i=l 



•occ), 



(6.4) 



where H KS is the Kohn-Sham Hamiltonian and p is the 

time-dependent density p(r,f) = E"" 1 i/ff(r,f)iA i (r,0- The 
solution of this equation relies on very simple sparse- 
matrix-vector multiplications and on a numerical imple- 
mentation of the unitary time evolution operator (Ya- 
bana and Bertsch, 1996, 1999a, 1999b; Bertsch, Iwata, 
et al, 2000; Bertsch, Rubio, and Yabana, 2000). Hence, 
for example, in the case of the TDLDA approximation 
for a cluster, the solution of the matrix equation (5.15) 
with a given kernel f xc is equivalent to propagating the 
equation above for some femtoseconds (the number de- 
pending on the accuracy in energy for the spectrum; Ya- 
bana and Bertsch, 1996, 1999a) with a given (now time- 
dependent) V xc . The real-time method has two major 
advantages: it is nonperturbative and therefore allows 
nonlinear effects of large fields to be calculated with the 
same effort, and it uses the same energy functional for 
the dynamical calculation as for the static calculation 
used to prepare the ground state, i.e., the potential and 
its density variation are automatically consistent, with- 
out the need to calculate the kernel explicitly. 

In principle, the Bethe-Salpeter equation can also be 
put into the same form, since the time-dependent den- 
sity can be obtained as the diagonal of the time- 
dependent Green's function (Kadanoff and Baym, 
1962). In that case, the time evolution equation for the 
Green's function must be solved. This has already been 
proposed by Kwong and Bonitz (2000) and applied to 
the model case of plasma oscillations, including damping 
effects on a correlated electron gas. The technique is 
similar in nature to solving the time-dependent Hartree- 
Fock equation (Langhoff et al, 1972); in principle, using 
the proper self-energy as in the GW approximation, the 
time evolution propagation is equivalent to solving the 
Bethe-Salpeter equation. An application to real systems 
has never been tried, although the idea seems appealing, 
since instead of solving a four-point equation one need 
only perform a sequence of (two-point) matrix-matrix 
multiplications, giving the action of 2 on G. However, 



since G is nonlocal, and since the time dependence of 2 
(and hence memory effects) cannot be neglected, it is 
not obvious whether such an application would be fea- 
sible, not least because of storage requirements. 



B. The limit of isolated electrons 

Before discussing results, it is useful to make a few 
remarks using a model system and/or simplified equa- 
tions. To this end, working with exchange only already 
allows one to understand many things. Note that, at least 
for finite systems, exchange only should be a much bet- 
ter approximation to reality for the two-particle excited 
states than for, say, photoemission, since screening is less 
important for a neutral than for a charged excitation. 

This exchange-only approximation has the advantage 
that one can work with relatively simple (explicitly 
known) exchange operators, and with the DFT exact- 
exchange potential (Gross et al., 1996). This latter po- 
tential (assuming from now on all wave functions to be 
real for simplicity) is given by 

V x (r) = 2lZ 2 J dr'j dr x J dr 2 ^i) 

M*i)<M* 2 ) , , , <Ac(r')<A„(r')\ 
x 2, — ltt^ti <Mr 2 ) — ; — ; — 



l r l- r 2l 



(6.5) 



where ^ (r,r') is the static independent particle re- 
sponse function. To simplify the following consider- 
ations further, let us suppose that the system has only 
one electron; V x is then (Taut, 1992) 
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The "quasiparticle correction" (ijj v \% x \ip v ) — (^ v \V x \ip v ) 
is hence zero for the occupied state, and for the lowest 
unoccupied state it becomes 
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Now, when considering only one occupied and the low- 
est empty state (single-pole approximation), the effect 
of the Hartree kernel and the electron-hole interaction 
is a change of the optical gap by the matrix elements of 
the Coulomb interaction v and the electron-hole inter- 
action (in this case also unscreened) taken with respect 
to the pair (v,c). The latter contribution is 

l^ c ( r)] 2 l^(r')l 2 

1 7] . ( 6 - 8 ) 

r-r 



i e—h . 



dr dr'- 



which exactly cancels the second term —{tf/ c \V x \ifr c ) of 
the quasiparticle correction [Eq. (6.7)], whereas the 
electron-hole exchange interaction (the Hartree kernel) 
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cancels the first term of Eq. (6.7), (if/ c \% x \if/ c ). On the 
other hand, in the TDDFT approach, the total correc- 
tion is simply the matrix element 

^ TDDFT _ J dr j dx' ljj*{x)lP v {x) 

which, using / A .(r,r') = 8V X I 8p(x')= — l/|r— r'| with the 
above V x [Eq. (6.5)], yields exactly the same result. 
Hence in both approaches the net correction to the ei- 
genvalue gap is zero, as it should be. This is of course 
due to the fact that in this simple model the valence 
exchange term is simply the self-interaction correction, 
and it completely cancels the Hartree term. Neverthe- 
less, it is instructive to study this case, since it illustrates 
how the matrix elements of f x simulate the sum of self- 
energy corrections and the electron-hole interaction. 
Moreover, even in more realistic systems, such as mol- 
ecules or clusters, the self-interaction contributions can 
be very strong. In particular, for other approximate ker- 
nels that do not treat exchange exactly, this cancellation 
can be incomplete and the one-electron limit may be 
wrong. 

If one admits more conduction states, it turns out 
again that in this exchange-only TDDFT all off-diagonal 
elements are zero, since the total kernel v+f x is zero. 
This is not true for the Bethe-Salpeter approach, where 
the contributions to the kernel coming from the Hartree 
part and the Fock part do not cancel, since they connect 
different indices. Therefore, in order to see whether the 
one-electron limit is correct, one should check the con- 
tribution of the off-diagonal elements. This is most eas- 
ily done in perturbation theory. Going beyond the diag- 
onal (first-order) approximation of the Bethe-Salpeter 
equation yields a second-order correction to the excita- 
tion energy which happens to be 

^(2) V \~{'Pc\^ x \>Pc')-{'Pc\Vh\>Pc')\ 2 in . 

E{'= 2j —qp qf: , (6.10) 

where the first term in the numerator is the matrix ele- 
ment between particle-hole pairs (v,c) and (v,c r ) of the 
electron-hole exchange, and the second term in the nu- 
merator is exactly equal to the matrix element of the 
direct electron-hole interaction. In the denominator, one 
has quasiparticle eigenvalues of empty states. This sug- 
gests a deviation from the correct one-electron result, by 
a shift — l-Ej^l to lower excitation energies. However, 
there are further corrections, due to the fact that the 
above formulas have been obtained (as is usually the 
case in a GW calculation) using Kohn-Sham wave func- 
tions, i.e., performing the GW calculation in first-order 
perturbation theory. One has in fact to calculate the cor- 
rection obtained by going to second order in the quasi- 
particle eigenvalues and taking into account the corre- 
sponding correction to the wave functions in the 
evaluation of the matrix elements of the exciton Hamil- 
tonian. It turns out that the change in eigenvalues adds 



another — |-E[ 2) | to the excitation energy, whereas the 
update of the wave functions shifts the optical gap by 
+ 2|_E^ 2 '|. Moreover, in this model system the resonant- 
antiresonant coupling terms are exactly vanishing. The 
final result is hence again correct. 

This suggests the following observation concerning 
TDDFT in the exact-exchange approximation and 
Bethe-Salpeter calculations for very small systems: It is 
reasonable to neglect off-diagonal elements in the 
TDDFT equation (single-pole approximation), as well 
as in the Bethe-Salpeter method, if Kohn-Sham wave 
functions are used and GW corrections are calculated to 
first order. However, if off-diagonal elements are in- 
cluded in the Bethe-Salpeter equation, it may be neces- 
sary to go beyond first-order perturbation theory in the 
GW calculation. On the other hand, resonant- 
antiresonant coupling terms can be neglected for those 
systems. 

One can also examine in detail whether the predic- 
tions made on the basis of these model calculations are 
valid. This can be deduced from the results of calcula- 
tions on SiH 4 done by Rohlfing and Louie (2000), in 
which Bethe-Salpeter results are compared to quantum 
Monte Carlo results (Grossman etal., 2001) and to ex- 
periment. First, Table III of Rohlfing and Louie (2000) 
confirms that the resonant-antiresonant coupling is in- 
deed completely negligible. Second, the same table con- 
firms that the inclusion of off-diagonal elements of the 
resonant electron-hole Hamiltonian lowers the excita- 
tion energy (by 0.41 eV for the first singlet transition). 
Then, Table II shows an increase of the transition energy 
by 0.67 eV due to the off-diagonal elements of X 
— V xc . This is in agreement with the above qualitative 
predictions concerning the cancellation of both effects. 
It suggests, moreover, as also pointed out by Rohlfing 
and Louie (2000), that most of what is going on can be 
understood on the basis of Hartree-Fock, self- 
interaction corrections and long-range potentials. In 
fact, the (not self-consistent) time-dependent Hartree- 
Fock results presented in the second column of Table II 
are clearly in acceptable agreement with the correspond- 
ing Bethe-Salpeter results. 

On the other hand, these findings should not be ex- 
trapolated to larger systems: in a solid, in fact, X pre- 
serves the crystal symmetry and can therefore not mix 
functions of different k points. The exciton Hamil- 
tonian, on the contrary, does strongly mix those transi- 
tions which are close in energy and which often come 
from the same band but different k points. Normally, 
one will find that quasiparticle wave functions are close 
to the Kohn-Sham ones, as has been pointed out by Hy- 
bertsen and Louie (1986), whereas a strong excitonic ef- 
fect that is entirely due to the mixing of transitions at 
different k points drastically alters the absorption spec- 
trum. Note that in a solid first-order perturbation theory 
with respect to the electron-hole interaction yields a 
vanishing shift of transition energies, and in fact the 
modifications in the spectra (even those appearing as 
shifts of peak positions) are essentially due to the mixing 
of transition matrix elements. As pointed out in Sec. 
IV.B, the joint density of states remains virtually un- 
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changed (Rohlfing and Louie, 1998b; Benedict and Shir- 
ley, 1999). Of course, this question about wave functions 
does not arise in the TDDFT approach: since the latter 
never pass through electron addition and removal ener- 
gies, only the Kohn-Sham wave functions resulting from 
the ground-state calculation are needed at all steps. 

To go further, one can now extend the model and con- 
sider the case of more than one valence orbital, but in 
the approximation that these orbitals are nonoverlap- 
ping. Again, one can start from the exact-exchange po- 
tential. Because the orbitals are nonoverlapping, one has 
v=v' in Eq. (6.5), which yields 



and 



X- 



\<l> v {*l)\ 2 *l>v>{*2)<l>c{*2)*l>c{*')ll>v{*') 



l r l" r 2l 

-1, 



X*o (*',*)• 



(6.11) 



Now, since the valence orbitals are nonoverlapping, 
Xo(r,r') is of block form in (r,r'), each of these blocks 
being given by a region R/ which contains both r and r' 
and having contributions from one valence orbital V[ 
only. We call these contributions xo v ,( r > r ')- The inverse 
of ^ is then given by the inverse of each block, and the 
above formula yields 



V x (r)=-J, f dt'dtrdt z 

V j ,Vj> J 



X- 



l^,(*i)| 2 

— : 7rXo Vl (r 2 y)Xo v >',*)• (6.12) 

l~ r 2\ ' 1 

The integral in r' has nonvanishing contributions of the 
form S(r 2 — r) when 1 = 1' and r is in Re, and one gets 



V x (r)=- dr lT — -j- 
IR. I r l~ r l 



(6.13) 



The integration is limited to the region R r where r is 
situated. The total quasiparticle correction stemming 
from both occupied and empty states is hence 

- ( f v \% x \ tfr v ) - ( &| V x \ i{, c ) + ( i/f v \V x \ ifr v ), 

„op / . iv i / \ i f j f i Pc(*)p(*i) 
A eP = (iA c S v <A C )+ dr\ dr x , , 

J jR r l*l _ *l 

f f ^,(r)^,(*)<M*')^(*') 

+ 2j " r " r i 



r— r 



f * f j P«( r )p( r i) 
dr dt] — j 

J JR T 1*1-* 



(6.14) 



In the diagonal approximation to the Bethe-Salpeter 
equation, the electron-hole exchange and the electron- 
hole interaction give, respectively, the contributions 

Ar / '= + J drj drV c (*)^(*)^(*')^(*')/|*-*'l 



dr\ dr'|^(r)| 2 |^„(r')| 2 /|r-r' 



as in the simpler model above. Adding these terms to 
A G/> , one finds that the total correction of the optical 
gap to the bare potential (T+V ) gap is 



AE° pt ~- 



dr dr x 



Pc(*)[p(*l)-P«(*i)] 



( the Coulomb interaction of the electron with the valence 
charge density after excitation), 



[ dr [ dr' 2 

J ■* v'*v 



r— r 



( the exchange interaction of the electron with the valence 
charge density after excitation), 

Pt,(*)p(*l) 



f i f , PoWpO 

- \ dr \ dri — j , 

J J l*i-*l 



(which subtracts the artificial Coulomb interaction of the 
now missing valence electron with valence charge density) 



dr J dr'2 



ip v {r)ip v ,{r) i/f v (r')i/f v (r) 



r— r 



( which subtracts the artificial exchange interaction of the 
now missing valence electron with the valence charge 
density, and, for the nonoverlapping orbitals, exactly can- 
cels the Coulomb self-interaction term). 

This result is very intuitive. In order to get the TDDFT 
result, one could now simply take the form (6.13) of V x 
derived above for the case of nonoverlapping orbitals, 
and obtain a kernel 



/*(*,*')= 



SV x (r) 
Sp(r') 



J R„R„ 



(6.15) 



Hence in the diagonal approximation the optical gap 
would turn out to be identical to the Kohn-Sham eigen- 
value gap, as in the case of only one electron, since f x 
cancels the Hartree part of the kernel. This is however 
not due to a failure of the exact-exchange approach; in 
fact, the kernel (6.15) has been derived from V x which 
was approximated for the case of the model of nonover- 
lapping valence states. Instead, one should first perform 
the functional derivative and then do the approximation. 
In that case, one gets additional terms, which should 
finally yield the correct result for the exchange-only 
case, as has been shown by Gonze and Scheffler (1999). 

It is also interesting to note that in the case of non- 
overlapping valence orbitals, the PGG kernel in the di- 
agonal approximation yields the same (wrong) result as 
Eq. (6.15). 

In summary, one should note that the exchange- 
correlation kernel has to simulate the result of self- 
energy corrections to the exact Kohn-Sham eigenvalues, 
plus the direct electron-hole interaction. In the one- 
occupied-level limit cancellations occur (for example be- 
tween the electron-hole interaction and the conduction- 
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state matrix element of V xc ) which might be used as an 
input for the design of kernels in systems with more 
electrons. From the case of several well-localized occu- 
pied states, one can see the importance of a kernel that 
is derived from a V x that well reproduces not only the 
ground-state density, but also density variations. 



C. Relaxation and correlation 

It is also worthwhile to analyze and compare the ef- 
fects of relaxation and correlation. For the N^N+ 1 
electrons excitation, one has of course a strong relax- 
ation (Hartree-like for localized states, and through cor- 
relation for delocalized states). Then, when going over 
from an N^N+ 1 to the N^N* excitation, one has to 
subtract the spurious electron-hole interaction and ex- 
change. Now there is also relaxation due to the presence 
of the hole. One can suppose that the relaxation due to 
the electron and the hole cancel to first order, so the 
total N^N* problem should have much less relaxation 
than the electron problem. 

The scenario describing exchange and correlation ef- 
fects in the one- and two-quasiparticle excitation spectra 
of many-body systems is hence the following: the screen- 
ing that explicitly appears in the GW self-energy comes 
from the electronic relaxation (classical or through cor- 
relation, as discussed in the introduction) due to the ad- 
ditional electron. The diagonalization of the GW Hamil- 
tonian takes into account the fact that the G W potential 
is different from the initial potential, which can be either 
Hartree, Hartree-Fock, V xc , or unknown if one does 
GW from scratch. Then, the diagonalization of the 
Bethe-Salpeter equation creates a correlated electron- 
hole wave function. Before the diagonalization, the 
electron-hole wave function is &{r e r h ) = ip*{r h )<p c {r e ) , 
which is uncorrelated. After the diagonalization, <t> has 
the correlated form <i> k (r e r h ) = Z vc A v k c i/,*(r h )if, c (r e ). In 
a molecule or cluster one can assume that one would get 
some reasonable result even by taking into account only 
diagonal elements; this is an uncorrelated, but interact- 
ing, electron-hole pair. It might also happen that only 
one occupied state is contributing, and the resulting 
wave function <^> x (r e r h ) = if/*(r h )[J, c A{i// c (r e )] would 
still be uncorrelated; only the electron has relaxed, i.e., 
adjusted itself to the fact that the hole is present. These 
mixing effects on the wave functions can in fact be easily 
visualized in a cluster (see, for example, Albrecht et al, 
1998b). In a solid, on the other hand, taking just one 
valence Bloch state would mean that one does not in- 
clude more than one k point, i.e., one would not get an 
excitonic effect. As in the case of the one-particle exci- 
tation, Hartree relaxation alone is not enough to de- 
scribe the solid: one must consider correlation. 

On the other hand, in TDDFT (and in the time- 
dependent Hartree approach) the off-diagonal elements 
immediately yield the correlated electron-hole pair, and 
the same discussion as above holds concerning the dif- 
ference between an extended crystal and a finite system. 



D. Comparison of some applications 

In real applications, of course, details of the system 
already at the independent-particle level are often the 
main ingredient for a first explanation of spectra. For 
example, bandstructure effects turn out to be decisive in 
determining not only the existence of interband transi- 
tions as observed in optical spectroscopy, but also the 
properties of plasmonic excitations in simple and noble 
metals. Indeed, ab initio calculations of the dynamical 
response have properly described the plasmon disper- 
sion in a variety of metals (Quong and Eguiluz, 1993; 
Aryasetiawan and Karlsson, 1994; Maddocks et al, 
1994a, 1994b; Fleszar et al, 1997; Ku and Eguiluz, 1999; 
Cazalilla et al., 2000). Since, however, the independent- 
particle approach often fails to yield quantitative, or 
even qualitative, agreement with experiment, it is useful 
to compare the performance of various approaches ap- 
plied to these realistic cases. 

In the previous section, we have shown or cited sev- 
eral examples for successful TDDFT calculations of the 
absorption spectrum of atoms and clusters. Far fewer 
examples exist for TDDFT calculations of the absorp- 
tion spectrum of bulk materials. The influence of an 
adiabatic LDA kernel f xc on the absorption spectrum of 
silicon has been calculated by Gavrilenko and Bechstedt 
(1996), who showed that its effect is negligible. These 
results have been confirmed by Koostra et al. (2000, al- 
though no comparison to the RPA result is given in that 
paper), and can also be seen in Fig. 10. In fact, since the 
independent-particle response function \o g° es to zero 
as q 2 , any kernel that does not have a II q 2 divergence 
does not contribute to the head (i.e., the G = 0, G'=0 
element) of the matrix f xc Xo and can at the best yield 
visible effects through off-diagonal elements, i.e., medi- 
ated by local field effects. Therefore it is often heard 
that "TDLDA works in clusters, but not in solids." We 
use this statement as a guideline for the following dis- 
cussions and to provide some explanations. In order to 
allow for a meaningful comparison, it is useful to con- 
centrate first on the effect of the part of the kernel which 
is common to the TDDFT and the Bethe-Salpeter ap- 
proaches, namely, the v (or v) part. In this way one can 
single out the rest and determine to what extent TDDFT 
and the Bethe-Salpeter approaches, respectively, actu- 
ally improve on P — the only quantity that is treated dif- 
ferently in the two approaches. 

1. Common ingredient: the bare Coulomb interaction 

As pointed out in Sec. II. B, the G=0 element of v 
determines the macroscopic screening l—v -Poo > which 
accounts for the essential difference between absorption 
and EELS spectra of solids. The rest of the Coulomb 
interaction, i.e., v, gives rise to what is called "local field 
effects" in the language of solid-state physics. This 
means that it reflects the fact that an inhomogeneous 
system will exhibit an electronic response that is position 
dependent (and not only distance dependent). It is intu- 
itively clear that such an effect will be the stronger the 
more a system is inhomogeneous. It is also obvious that 



Rev. Mod. Phys., Vol. 74, No. 2, April 2002 



644 



Onida, Reining, and Rubio: Density-functional vs many-body 



an atom, a molecule, or a cluster of whatever size is by 
itself a strong inhomogeneity in the empty space, and 
that the electronic response must depend strongly on the 
distance and on the position of the perturbation and the 
probe. 46 Therefore one can expect that the effect of v 
alone will already be a strong modification of the spectra 
in a cluster, whereas in a solid it will depend on details of 
the charge density of the latter, with a tendency to have 
very weak effects in, say, nearly free electron metals (ho- 
mogeneous systems) like simple metals. This is in fact 
the case. In Fig. 8 we have already shown the results for 
three different clusters, and discussed the fact that the 
Coulomb term has a very strong effect and removes 
most of the discrepancy between the LDA curves and 
experiment. 

When moving to solids, a good example is the loss 
function of bulk silicon, in which the localization of the 
charge on the bonds gives rise to a non-negligible inho- 
mogeneity. The results of local field effects on the plas- 
mon resonance of silicon have already been discussed, 
for example, by Louie et al. (1975), who found that the 
local field effects are not as pronounced as in the case of 
the cluster, but are sizable. Moreover, as in the case of 
the cluster, agreement with experiment is improved 
when v is included; in particular, the height of the plas- 
mon peak is considerably lowered. 

When local field effects are included, good-quality 
loss spectra are also obtained in transition-metal com- 
pounds like Ni and NiO (Aryasetiawan, 1994). The good 
quality of RPA loss function results and the importance 
of local field effects is illustrated for rutile Ti0 2 in Fig. 
11 (Vast et al, 2002). One can see that local field effects 
are particularly important at higher energies, in the re- 
gion of transitions from the Ti 3p semicore level, yield- 
ing good agreement with experiment. This almost seems 
to suggest that one could, in a first approach, always 
neglect the cumbersome exchange-correlation effects. 
This is of course not true. It is instructive to look, not 
only at the loss spectrum, but also at the absorption 
spectrum of bulk silicon [Im(e M )], shown in Fig. 5. The 
dot-dashed curve (RPA result) is in fact in poor agree- 
ment with experiment. Another striking example is the 
absorption spectrum of solid argon, shown in Fig. 12. 
The strong absorption peaks in the experiment (Saile, 
1976; Saile et al., 1976) on the low-energy side, which are 
known to be of excitonic nature, are completely missing 
in both the RPA and the GW -RPA results (Olevano and 
Reining, 2001b). Note that local field effects are in- 
cluded in these calculations, but they cannot remove the 
discrepancies. The effect of v on Im(e M ), in its more or 
less pronounced form, is in fact essentially to shift oscil- 
lator strength to higher energy. This is due to the fact 
that v is positive, being the Coulomb interaction be- 
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46 The essential difference between even a very big cluster and 
an infinite solid is the presence of a surface between the cluster 
and vacuum, implying boundary conditions for the electric 
field. 
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FIG. 11. Integrated loss function at q^OA A -1 for rutile Ti0 2 
(Vast et al., 2002): solid line, experiment; dashed line, RPA 
results. The dot-dashed line in the inset shows the results with- 
out local field effects. 

tween charge-density waves (Hanke and Sham, 1975). It 
never creates or even significantly enhances structures 
on the low-energy side. 

Hence v, or local field effects, are generally not suffi- 
cient to obtain quantitative (for e _1 ) or even qualitative 
(for e M ) agreement between theory and absorption ex- 
periments, and one has to add the exchange-correlation 
effects. In other words, one has to use the interacting P 
(instead of P = P IQP ) in Eqs. (2.10) and (2.11). With re- 
spect to the RPA-plus-local-field result, these effects 
should cause a variety of modifications in the spectra, if 
they are supposed to restore agreement with experi- 
ment. Essentially, from the above results one can see 
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FIG. 12. Optical spectra for argon: •, experiment (Saile, 
1976); dashed line, RPA; dot-dashed line, TDLDA; solid line, 
GW-RPA. From Olevano and Reining, 2001b. 
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FIG. 13. Dynamical structure factor S(q,&>) 
= — 2ft fdrdr'e~ ,q(r ~ r ' Im x(r,t' ,co) of aluminum, calculations 
from Fleszar et al. (1995). Open circles are experimental data 
from Platzman et al. (1992). In the left panel, they are com- 
pared with the Kohn-Sham independent particle response (full 
line labeled S KS , obtained taking x = Xo in the above expres- 
sion for S) (the double -bump structure was already obtained 
by Maddocks et al, 1994a, 1994b). The dotted line, labeled 
Sundhardi is the corresponding noninteracting response of jel- 
lium with r s = 2.07. In the right panel, the same experimental 
data are compared with S calculated using x from Eqs. (5.11) 
and (5.12). Full line, TDLDA (f xc ) from Eq. (5.13); dashed 
line, RPA (f xc = 0). 

that a shift is needed, which can be towards either lower 
energies (as in the case of the loss spectra of both clus- 
ters and solids) or higher energies (as in the case of the 
absorption spectrum of bulk silicon), together with a re- 
distribution of oscillator strength towards lower ener- 
gies, which may enhance existing peaks and/or create 
new ones, as in the case of solid argon (see Fig. 12). 

2. Exchange and correlation effects 

Concerning these exchange-correlation effects, the 
Bethe-Salpeter and the TDDFT approaches act in a 
completely different way. In fact, the Bethe-Salpeter ap- 
proach can be seen as a two-step method, where first a 
GW calculation yields a strong shift of the whole spec- 
trum towards higher energies, and the subsequent inclu- 
sion of the electron-hole interaction redistributes oscil- 
lator strength towards the low-energy peaks, and 
eventually shifts peaks, or even creates peaks in the qua- 
siparticle gap. It is clear that there must be a partial 
cancellation of the self-energy corrections and the 
electron-hole interaction. In the small cluster discussed 
above, and also for bulk plasmons, where the RPA plus 
local fields calculation already gives very good results, 
this cancellation is almost perfect. TDDFT, on the other 
hand, should describe directly the complete electron- 
hole excitation, but in an effective way, which may be 
more difficult to understand [see also the discussion in 
Appendix A, after Eq. (A5)]. In this case, too, partial 
cancellations occur. Figure 13 shows an example of these 
cancellations. Fleszar et al. (1995) have calculated the 
dynamical structure factor of bulk aluminum within 
TDDFT. The independent-particle response ^ closely 
simulates the measured spectrum, and even the double- 



hump structure is correctly predicted. However, when a 
full RPA calculation, including the Coulomb kernel, is 
performed, the result worsens. Exchange-correlation ef- 
fects should hence completely cancel the Coulomb re- 
pulsion in that case (i.e., a kernel f xc = — l/|r— r'| would 
do the job). Depending on the situation and the system, 
different aspects of f xc will become important. For ex- 
ample, Sturm and Gusarov (2000) have shown that in 
order to describe the aluminum inelastic x-ray scattering 
cross section for large momentum transfer and large w, 
dynamical correlation effects in f xc are more important 
than band-structure effects (interband transitions are in- 
sufficient to explain the observed magnitude of the high- 
energy tail of the spectrum, while multiple particle exci- 
tations give the correct order of magnitude). 

a. Some applications to finite systems 

For clusters, the majority of DFT-based calculations 
today rely on the use of the static LDA kernel, since the 
optical spectra turn out to be in very good agreement 
with experiment even at the TDLDA level (Casida, 
1995; Rubio et al., 1996; Yabana and Bertsch, 1996; Va- 
siliev et al, 1999, 2002). Although, as pointed out above, 
the improvement with respect to a static LDA calcula- 
tion essentially comes from the Coulomb part of the ker- 
nel, and not from f xc , these findings have led to an im- 
portant breakthrough in the calculation of excitation 
spectra of finite systems and make the TDDFT tech- 
nique very popular. One example is the benzene mol- 
ecule (Fig. 14). It is clear that apart from small differ- 
ences due to the resolution of vibrational modes in the 
experimental spectrum, the overall response is very well 
described by the TDLDA approach. 47 Satisfactory 
agreement is also obtained for other molecular clusters 
as shown in Fig. 8. This figure however, also shows the 
need for improvements over the simple TDLDA. 

A good and widely studied example is the silane mol- 
ecule SiH 4 , for which quantum Monte Carlo (Grossman 
et al, 2001), Bethe-Salpeter (Rohlfing and Louie, 1998a, 
2000; Grossman et al, 2001), and TDDFT results (Ogiit 
et al, 1997; Marques et al, 2001) are available. In par- 
ticular, TDLDA or gradient corrected functionals pre- 
dict that bound excitations will be resonances. This arti- 
fact is clearly corrected by the exact-exchange potential 
and a method proposed by van Leeuwen and Baerends 
(1994) to impose the correct asymptotic behavior on the 
exchange potential (Marques et al, 2001). Within this 
scheme the TDDFT results for SiH 4 have the same ac- 



For the chiral C 76 fullerene, TDLDA also provides a con- 
sistent description of the optical absorption spectrum, the cir- 
cular dichroism spectrum, and the optical rotatory power (ex- 
cept for an overall shift of the total spectrum; Yabana and 
Bertsch, 1999a, 1999b). Similarly, the optical spectrum can be 
used to elucidate the ground-state structure of the smallest 
possible carbon cage C 2 o (Castro et al, 2001), which has been 
elusive to ground-state calculations in both quantum Monte 
Carlo (Grossmann et al, 1995) and DFT (Jones and Seifert, 
1997). 
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FIG. 14. Optical absorption of the benzene molecule, in units 
of eV~ 1 . The TDLDA calculation is from Yabana and Bertsch 
(1999a), whereas the experiment is from Koch and Otto 
(1972). 

curacy as the GW or quantum Monte Carlo 
calculations. 48 This fact is illustrated in Fig. 15. 

TDLDA gives good results for Na 4 as well as for small 
and large clusters of simple and noble metal atoms (Ru- 
bio and Serra, 1993; Rubio et al, 1996, 1997; Yabana and 
Bertsch, 1996; Serra and Rubio, 1997; Vasiliev et al, 
1999, 2002) and fullerenes (Yabana and Berstch, 1999a, 
1999b; Castro et al, 2001) but works less well in small 
hydrogenated silicon clusters (see Fig. 15) 49 This prob- 
lem might be due to the presence of more localized 
bonds. The exact-exchange calculation does help in this 
case and keeps nearly the same accuracy as in the so- 
dium cluster. Both calculations are in slightly better 
agreement with experiment than the Bethe-Salpeter cal- 
culation (Onida et al, 1995). Note that the two peaks at 
about 3 eV in both the exact-exchange and the Bethe- 
Salpeter calculations are very similar. However, in both 
clusters we see an overall tendency of the TDDFT cal- 
culations with approximate kernels to give larger excita- 
tion energies (peak positions blueshifted with respect to 
experiment). By looking carefully at the different contri- 
butions coming from the correct asymptotic behavior of 
the potential, as well as at variations of the f xc kernel, 
Marques et al. (2001) show that there is an inherent 



48 Very good results are obtained in the three approaches — 
quantum Monte Carlo, Bethe-Salpeter and TDDFT with exact 
exchange. In particular, they overestimate the energy of the 
first singlet excitation by 3%, 4.5%, and 1.5%, respectively, 
whereas the ionization potential (exactly reproduced in quan- 
tum Monte Carlo) is only 4% overestimated by TDDFT with 
exact exchange. The quantum Monte Carlo results are equiva- 
lent to the best quantum-chemical (coupled-cluster and com- 
plete active space self-consistent-field) calculations (for details 
see Grossmann et al, 2001 and references therein). 

49 Here it is important to remark that for the optical gaps, the 
calculated values within TDLDA for large-size hydrogenated 
silicon clusters (Ogiit et al, 1997; Vasiliev et al, 2001) agree to 
within =10% with experiment. 
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FIG. 15. Calculation of the strength function for Na 4 and SiH 4 
within TDDFT using different kernels (Marques et al, 2001): 
solid line, TDLDA; dotted line, exact-exchange. Gray solid 
line, experimental results; dashed line, the Bethe-Salpeter 
equation results of Onida et al. (1995) for Na 4 . Vertical lines 
are the lowest excitation energies obtained by Grossman et al. 
(2001) for SiH 4 : solid vertical line, quantum Monte Carlo cal- 
culations; dashed vertical line, Bethe-Salpeter equation. In the 
inset: solid line, self-interaction correction; dotted line, van 
Leeuwen and Baerends (LB94) prescription. 

source of error in the approximated kernels coming 
from neglect of dynamical (correlation) effects. These 
tend to shift the spectrum to lower energies and are also 
responsible for improving the excitation energies of 
Rydberg-like states (where exchange effects are small 
due to the nearly zero overlap between the Rydberg un- 
occupied state and the other valence states; thus the usu- 
ally weak polarization effects are going to play an im- 
portant role here). 

For these small molecules, there are thus various 
methods to get good transition energies: quantum 
Monte Carlo, the Bethe-Salpeter equation, and TDDFT 
yield satisfactory results. One might of course wonder 
whether the calculation of spectra involving empty (ex- 
tended) states is not questionable in a supercell ap- 
proach. In fact, high-energy, extended states do repre- 
sent a problem if they are really needed as the physical 
final (one-electron) state (e.g., in absorption at higher 
frequency, or for a photoemitted electron, if one wants 
to go beyond the usual approximation of neglecting the 
final-state effects). Most often, however, these states are 
used only as intermediate states in summations, i.e., for 
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a closure relation, which is a mathematical property de- 
riving from the completeness of the basis set. Since the 
validity of the closure relation does not depend on the 
chosen boundary conditions (even if the actual shape of 
the individual states does), this does not present a prob- 
lem. 

Concerning the Bethe-Salpeter approach for SiH 4 , it 
should be noted that (i) the results depend on the 
choices made for the different steps of this rather com- 
plicated technique, e.g., the diagonalization of the GW 
Hamiltonian discussed above; and (ii) only transition en- 
ergies, not line shapes, are shown. However, in TDDFT 
the results depend of course on the choice made for the 
kernel, as discussed above. 

b. Some applications to extended systems 

In solids one should look at both absorption and 
electron-energy-loss spectra. As discussed in the previ- 
ous sections, TDLDA fails to describe the absorption 
spectra of solids. In fact, at least in the case of silicon, 
the failure of the TDLDA f xc to reproduce the correct 
long-range behavior of the kernel explains its bad per- 
formance, as shown above in Fig. 10. On the other hand, 
Olevano and Reining (2001a) have calculated the 
electron-energy-loss spectrum of silicon including local 
field effects, GW corrections, and excitonic effects, for 
vanishing momentum transfer. The inclusion of excitonic 
effects improves the results with respect to both the 
RPA and, more drastically, the GW calculations (see 
Fig. 16): GW and excitonic corrections cancel to a large 
extent. In contrast to the case of small molecules and to 
absorption spectra in bulk, here the correct result could 
only be obtained by taking into account the coupling 
between transitions at positive and negative frequencies 
[blocks K in Eq. (B20)]. Figure 16 also shows that 
TDLDA improves with respect to the RPA calculation 
(including local field effects), which means that even the 
TDLDA f xc has a visible effect on the electron-energy- 
loss spectrum. However, as in the case of clusters, local 
field effects are more important than f xc for the 
electron-energy loss spectra. This is a very general find- 
ing in TDLDA electron-energy-loss calculations (see the 
results of Waidmann et al, 2000 and references therein). 

For nonvanishing momentum transfer, one can com- 
pare the TDDFT results of Waidmann et al. (2000) on 
diamond to the Bethe-Salpeter results of Caliebe et al. 
(2000). The comparison is limited, however, to the low- 
energy region well below the plasmon. First of all, the 
work of Waidmann et al. (2000) shows the increasing im- 
portance of local field effects with increasing q. Local 
field effects start to be visible at q = l.QAT l (about 1.1 
YX). At g = 1.5A _1 (about 1.7 YX) they are already 
very strong. Second, TDLDA gives a reasonable de- 
scription of the low-energy part of the spectrum. It 
should be pointed out that the authors state that the 
LDA kernel f xc itself has only a negligible effect on this 
result. On the other hand, in the work of Caliebe et al. 
(2000) the Bethe-Salpeter method is applied to the cal- 
culation of the same spectra. Results of comparable 
quality are obtained in the end. A more detailed com- 
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FIG. 16. Energy-loss spectra of bulk silicon: •, experiment 
(Stiebling, 1978); dotted line, Bethe-Salpeter equation without 
resonant-antiresonant coupling; dash-dotted line, GW-RPA; 
dashed line, RPA; dash-double-dotted line, TDLDA; solid 
line, Bethe-Salpeter equation with coupling (Olevano and 
Reining, 2001a). 

parison is difficult, since the Bethe-Salpeter results are 
compared to those obtained when GW corrections are 
added to the Kohn-Sham eigenvalues in the 
independent-particle polarization. Moreover, as is often 
the case in the Bethe-Salpeter approach, the local field 
effects are treated as part of the electron-hole interac- 
tion kernel. In other words, part of the electron-hole 
effect shown is due to the local field effects. The net 
improvement due to the electron-hole interaction is 
hence very strong (although not systematic for all spec- 
tra), but it is difficult to estimate the total improvement 
due to self-energy plus electron-hole attraction effects, 
in comparison to RPA (including local field effects), or 
better, TDLDA, results. A detailed comparison of 
TDDFT and Bethe-Salpeter results for the q # case is 
in fact still missing in the literature. Bethe-Salpeter re- 
sults should be superior to TDLDA results in those 
parts of the loss spectra that are dominated by interband 
transitions, i.e., by the structures in the imaginary part of 
the dielectric function, which are not well reproduced in 
TDLDA. 

VII. CONCLUSIONS— FREQUENTLY ASKED QUESTIONS 
AND OPEN QUESTIONS 

Since the very beginning of physics and chemistry, the 
interaction of photons and electrons with matter has 
been a major topic of study. Today, most characteriza- 
tion tools as well as electro-optical devices are based on 
our understanding of these interactions. Technological 
applications are rapidly progressing, but many funda- 
mental questions concerning theoretical and numerical 
descriptions are still open. 
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This has motivated us to discuss in the present review 
the two most widely used techniques for describing elec- 
tronic excitations in finite and infinite systems, namely, 
the Green's-function approach to many-body perturba- 
tion theory calculations, and time-dependent density- 
functional theory. We have presented the foundations 
and approximations of the two approaches as well as the 
advantages and drawbacks of the methods. A special ef- 
fort has been made to provide a unified picture of the 
underlying equations, in order to facilitate the compari- 
son and to analyze the most frequently used approxima- 
tions. We have stressed the fact that the two approaches 
are somewhat complementary and that one should work 
on both sides in order to optimize the search for system- 
atic improvements. 

We have not presumed to give a complete view of all 
theoretical descriptions of spectroscopies, but rather 
looked at the field from a particular point of view, with 
the aim of motivating some possible future develop- 
ments. Some points have been treated in more detail 
than others; we have made our choices on the basis of 
the very lively discussions which are now going on in this 
field. As a conclusion, we summarize below a set of fre- 
quently asked questions, which we have tried to answer 
in the present paper. Many points are not new, and they 
will be obvious for a specialized reader, but we think 
that it is useful to put them in a common context. This 
should also be helpful to those who wish to enter the 
field. 

Other questions do not have an answer yet. They will 
be treated at the end of this section — as a motivation for 
future research. 

• The delta-self-consistent field approach, i.e., explicitly 
calculating total energy differences, is known to work 
well. Why can't we just always calculate total energy 
differences in DFT, instead of struggling with the mean- 
ing of Kohn-Sham eigenvalues! First, ASCF in DFT 
supposes that one can simulate the initial and the final 
state of the excitation by occupying selected one- 
particle orbitals. This excludes excitations that are not 
easily described in terms of isolated single -particle 
transitions (for example, the collective plasmon exci- 
tations). Second, whereas the choice of the N and N 
± 1 states can be straightforward in a small system, it 
is not necessarily well defined in the bulk; in fact, 
when Bloch states are chosen to describe the elec- 
trons, Hartree relaxation effects vanish and hence the 
main advantage of ASCF is lost. One should then in- 
clude dynamical correlation effects, which go beyond 
static DFT-ASCF. See discussions in Sees. IV.A.2 and 
VI.C. 

• Does TDDFT give electron removal or addition ener- 
gies! In Sec. V we described how TDDFT provides an 
exact framework for getting the excitation energies of 
an iV -particle system. The theory handles only neutral 
excitations, that is, excitations in which the number of 
particles does not change. Concerning N^N±1 pro- 
cesses, TDDFT is only supposed to have the correct 
ionization potential threshold, as this should be guar- 



anteed by the exact DFT exchange-correlation poten- 
tial (Almbladh and von Barth, 1985). 

• If quasiparticle energies are electron addition and re- 
moval energies, i.e., total energy differences, how can 
they be complex! In the Lehmann representation of 
the Green's function the excitation energies involved 
are always differences of total energies between the N 
and N±l systems; hence they are always real. The 
infinitely close-lying peaks of the spectral function in 
the Lehmann representation merge into broad struc- 
tures that can be identified as quasiparticle peaks. 
These structures have a finite width and can hence be 
described by a real part (peak position) and an imagi- 
nary part (width) of a quasiparticle's energy. See dis- 
cussion in Sec. IIIA. 

• Why can TDDFT yield absorption spectra, when DFT 
is a ground-state theory! DFT is a ground-state theory 
because it is based on a minimization of the total en- 
ergy. TDDFT, on the other hand, is derived from the 
extrema of the quantum-mechanical action, and hence 
describes the evolution of the system instead of its 
ground state. This difference is equivalent to that in 
classical mechanics, where the equilibrium position of 
a particle in space can be found by looking for the 
minimum of the potential, whereas its trajectory is de- 
rived from the extrema of the classical action. See dis- 
cussion in Sec. V. 

• Can one in principle get excitonic effects in TDDFT! 
Yes, excitonic effects are in principle exactly contained 
in the TDDFT equations, if the exact exchange- 
correlation potential V xc and kernel [f xc (r,r' ,a>)] are 
used. However, the practical implementations use 
simple functionals, which lack, for example, the 
proper spatial nonlocality. Therefore direct electron- 
hole interaction effects are only partially described in 
finite systems, and in general still out of reach of to- 
day's TDDFT calculations for solids. See Sec. V.C and 
Appendix C. 

• Is the fact that Kohn-Sham eigenvalue differences un- 
derestimate the gap a problem for the calculation of 
absorption spectra within TDDFT! No. The Kohn- 
Sham eigenvalue differences are not meant to repro- 
duce the (optical) gap. The optical gap is determined 
by the Kohn-Sham potential (eigenvalues) and by 
variations of the Kohn-Sham potential, through the 
kernel. This latter contribution changes the optical 
gap with respect to the Kohn-Sham eigenvalue differ- 
ence. Of course, results improve when better poten- 
tials (eigenvalues) are used — but "better" does not 
mean "close to quasiparticle energies." 

• Why does the RPA often yield good results for 
electron-energy loss, but not for absorption in solids! 
This can be traced back to the crucial role played by 
the long-range part of the Coulomb potential in the 
kernel: the exchange-correlation contribution f xc is in 
fact added to the Hartree term v in the case of the loss 
spectra [see Eqs. (2.3) and (2.11)], but to the differ- 
ence between the bare v and the long-range part 
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u(G = 0) in the case of absorption [see Eqs. (2.9) and 
(2.10)]. In the latter case, the only long-range term 
comes from f xc , which makes it so crucial. See Sees. 
II.B and VI.D. 

• Why does TDLDA work for clusters, but not for the 
absorption spectra of solids? First, the main difference 
between the results of TDLDA and those of Eq. (2.8) 
[the IP-RPA dielectric function] in clusters stems from 
the contribution v in the kernel, which is always exact. 
Second, for finite systems f T x ® LDA provides an effec- 
tive electron-hole attraction, which gives a significant 
contribution to the spectra. In infinite nonmetallic sys- 
tems, the LDA/^. only contributes through local field 
effects, since the head of the matrix f xc Xo vanishes as 
q^Q. See Sees. V.C, VI.D, and Appendix C. 

• What is the link between the Hartree potential, the 
electron-hole exchange interaction, and local field ef- 
fects in a solid! The variation of the Hartree potential 
with respect to the density yields a contribution v to 
the kernel of both the TDDFT and the Bethe- 
Salpeter equations. This contribution has to be taken 
with or without (v or v) its long-range (G = 0) contri- 
bution in the case of electron-energy loss or absorp- 
tion, respectively. Neglecting v in the screening equa- 
tions is equivalent to neglecting local field effects. The 
matrix elements of v and v in an electron-hole basis 
are dipole-dipole like, and one calls this contribution 
therefore the "electron-hole exchange" in the Bethe- 
Salpeter framework. See Sees. II.B and IV.B.2 for 
more details. 

• Is the electron-hole exchange interaction screened or 
unscreened? The electron-hole exchange interaction 
stems from the density variation of the Hartree poten- 
tial and is therefore exactly the unscreened Coulomb 
interaction (see Sees. II.B, II. C, and IV.B.2). Some- 
times model calculations can be found in the literature 
where this interaction is screened. This procedure is 
justified when (and only when) an extremely re- 
stricted space of transitions is explicitly mixed by the 
exciton; the neglected transitions are then implicitly 
put back into the calculation by an effective screening 
of the electron-hole exchange. 

• To describe two-particle excited states one resorts to the 
two-particle Bethe-Salpeter equation. What about 
three-, four-, or more-particle excitations? Do we need 
three- or more-particle analogies to the Bethe-Salpeter 
equation? No. In fact, even for the case of two-particle 
excited states we do not need the two-particle Bethe- 
Salpeter equation, since it would be formally sufficient 
to solve Hedin's equation for the vertex [Eq. (3.25)]. 
The situation is similar to the case of electron addition 
or removal energies, where the set of higher-order en- 
tangled Green's functions is decoupled by introducing 
the concept of the self-energy, which implicitly con- 
tains all higher-order electron interactions. Similarly, 
the complex equation for the polarization operator in 



terms of higher-order polarization operators can be 
formally solved by introducing the vertex function. 
See Sec. IVB. 

• TDDFT equations can be written as a quadratic equa- 
tion (6.3) or as a Dyson-like equation (5.11): are the 
two forms equivalent, and what are the reasons to 
choose one of them? Yes, the two forms are indeed 
equivalent. The TDDFT equation can be written for 
four-point response functions or for two-point (con- 
tracted) functions [see Eqs. (2.15) and (2.19)]. The 
two-point formulation (which does not exist for the 
Bethe-Salpeter equation) leads to the Dyson-like 
equations for two-index matrices [Eq. (5.11)], with the 
obvious advantage of a better scaling with the system 
size. On the other hand, the four-point formulation 
allows one to switch easily to transition space and to 
compare directly with the Bethe-Salpeter equation. 
The transition-space formulation can be written as an 
eigenvalue equation like Eq. (5.15) or as a quadratic 
equation for co like Eq. (6.3). It is advantageous for 
cases in which only a limited number of transitions 
contribute to the spectrum. See Sees. VB, VIA, and 
II.C. 

• Since most calculations use supercells or finite do- 
mains, can we trust the higher-energy eigenstates, which 
can "sense" the boundary conditions? If one is explic- 
itly interested in these states, that can actually be a big 
problem. However, most often (e.g., in the calculation 
of a Green's function) they are just virtual states that 
are summed over, and one does not care. See also Sec. 
VI.D.2.a. 

One important question now is actually which way 
one should go — is the Bethe-Salpeter equation or 
TDDFT more promising? Since both approaches are ex- 
act, the answer depends of course on whether — and 
when — the remaining open questions in the two ap- 
proaches can be solved. Some questions are common to 
both. An important point is whether pseudopotentials 
are always adequate for describing spectra with the pre- 
cision one asks for today. For certain applications, one 
should perhaps migrate towards all-electron schemes. 
Other questions are specific to the Bethe-Salpeter or the 
TDDFT approach. 

The Bethe-Salpeter approach offers a clear physical 
picture and straightforward possibilities for the analysis 
of results. It seems to work over a wide range of systems. 
The TDDFT approach, on the other hand, is appealing 
since it calculates things in a more direct way (without 
passing through electron addition and removal energies) 
and is, in principle, easier to use. 

Open problems in the Bethe-Salpeter approach in- 
clude, of course, the need for technical developments to 
overcome the bottlenecks linked to the four-point equa- 
tions. Moreover, just as in G W, not all the "ingredients" 
of the method are uniquely defined: how, for example, 
should vertex corrections and dynamical screening be 
included consistently? It is clear that the level of ap- 
proximation which should be used for each ingredient in 
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the various Green 's-function-based approaches is a deli- 
cate and nontrivial question. One should wonder why in 
the calculation of an electron-energy-loss spectrum it is 
a good approximation to use e _1 constructed from a P 
= Pjq P with DFT eigenvalues, and a much worse one 
when GW eigenvalues are used, whereas in GW calcu- 
lations, especially for large-gap systems, an update of 
the energies entering in e ~ 1 yields improvements for the 
resulting quasiparticle energies. At first sight, this seems 
contradictory. At this point we propose that one answer 
may be to put all choices on the basis of a consistent 
iteration scheme. This means that one starts with a GW 
calculation and subsequently iterates Hedin's equations. 
The iteration first yields P, which contains a vertex cor- 
rection. This explains why the electron-energy-loss spec- 
trum using P = Pjq P with GW eigenvalues yields poor 
results, but using P with GW eigenvalues and a vertex 
correction yields good results. Then one can go further 
and put this P into the equation for X. Now a second 
vertex is appearing, and it is has been shown that there 
are cancellations between the vertex in P and the ex- 
plicit one in X. It is known that it is better to consistently 
neglect the vertex in both W and X, which in practice 
means that a simple update of the energy denominator is 
in general the most consistent thing to do in a GW cal- 
culation, if one does not want to include both vertices. 
Note that one quickly runs into trouble when also up- 
dating the wave functions, because of the dynamical ef- 
fects, whereas it sometimes turns out to be necessary to 
update their spatial behavior. These are handwaving ar- 
guments, though essentially these — and the success of an 
approach in real calculations — justify the choices that 
are made in the iteration of the equations. More 
successes — and failures — of the Bethe-Salpeter ap- 
proach in its actual form are needed in order to sort out 
this question. 

In TDDFT, calculations are less cumbersome (if the 
two-point representation is chosen) than in the Bethe- 
Salpeter method at the moment, i.e., using the simple 
kernels available today, like that of the adiabatic LDA. 
Better exchange-correlation potentials and better ker- 
nels have to be found, especially if TDDFT is meant to 
become a method for the calculation of absorption spec- 
tra in solids. These kernels might well turn out to be 
very complicated, and as difficult to treat as the Bethe- 
Salpeter equation. There is some hope that, at least in 
certain cases, relatively simple solutions can be found 
(see Appendix C). The Bethe-Salpeter equation seems 
to provide a good starting point for the derivation of 
such effective kernels, since it contains the essential 
physics in a structure that is actually close to that of 
TDDFT 

In conclusion, it seems reasonable to suppose that 
progress will come from a common effort of people 
working in the fields of the Bethe-Salpeter equation and 
of TDDFT. 
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APPENDIX A: EXACT PROPERTIES OF f xc 

Although the development of time-dependent func- 
tional is still at a very early stage compared to that of 
static functionals, some of the known results for a homo- 
geneous gas can be generalized to the inhomogeneous 
case. For example, causality leads to Kramers-Kronig 
relations for the real and imaginary parts of f xc (r,r' 
and the fact that f xc (r,t;r' ,t') is a real-valued quantity 
implies that f xc (r;r' ,co)=f xc (r;r' ,— co)*. Besides that, 
the response functions satisfy the spatial symmetry rela- 
tions ^(r,r',a>) = ^(r',r,<t)), provided that the unper- 
turbed system has time-reversal symmetry (no magnetic 
field is applied or generated in the system). Thus one 
finds that f xc (r,r ' , to) =f xc (r ' ,r, w) . 

Further constraints on the potential V xc and kernel f xc 
can be deduced from the quantum-mechanical equation 
of motion applied to the position operator f: 
(d/dt(V(t)\r\y(t)) = i(V(t)\[H(t),i]\V(t))). The final 
and rigorous result is (Gross et dL, 1996) 

J drp(r,t)VV xc [p](r,t) = 0. (Al) 

This relation is also satisfied by the Hartree potential. 
This is nothing other than a mathematical statement of 
the physical fact that the exchange-correlation potential 
does not exert a net force on the system. One also ob- 
tains the corresponding "zero-torque" expression by 
considering the angle operator ip and using the rota- 
tional invariance of the Coulomb interaction, yielding 

J dip(r,f)rXVV«[p](l,f) = 0. (A2) 

Corresponding properties of the exact exchange- 
correlation kernel are obtained by evaluating the two 
previous expressions at the density p(r,f) = p (r) 
+ Sp(r,t) for arbitrary Sp(r,t). In frequency space one 
gets 
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rfrp (r)V/ xc { Po ](r,r>)=-V r ,L x Jp ](r'), (A3) 



rfrp (r)rxV/ xc [p ](r,r',a,)=-r'xV r ,y. YC [ Po ](r'). 

(A4) 

Concerning the electronic damping mechanism which is 
included in f xc , we should emphasize that the standard 
static approximations for f xc [where / xc (w):=/ lc (w = 0) 
is a real quantity] based on the homogeneous gas predict 
an infinite lifetime for plasmon excitations at small wave 
vectors. This artifact is inconsistent with the fact that the 
exact kernel has an imaginary part related to the multi- 
pair component of the response function Xmpil^) 
(Sturm and Gusarov, 2000), that is, 



Im[/* c (q,ft>)] = 



u(q) 2 Imx (q,w), 



where w p is the plasmon frequency. These higher-order 
interactions are beyond the usual description of screen- 
ing within the RPA or bubble approximation, and in- 
volve the partial summations of many complex diagrams 
appearing in the many-body response function. 

In the homogeneous electron gas some exact relations 
for the f xc kernel must hold. In particular one has the 
following relations: 

(i) the compressibility sum rule, lim q ^ f xc (q ,w = 0) 
= d 2 pe xt .(p)/dp 2 , where e xc (p) denotes the 
exchange-correlation energy per particle; 

(ii) the third frequency moment sum rule: 



0= 



,2/3 



dp 



2/3 



+ 6p 



1/3 



de xc {p)lp 
d p 



1/3 



(iii) the static and frequency-dependent short- 
wavelength behavior: \im (j ^ x f xc (q,co = 0)^c 
— (b/q 2 )[l— g(0)] [the correct values of the con- 
stants b and c can be found in Moroni et al. 
(1995) instead of the values used by Gross and 
Kohn (1985)] and lim^ x / xc .(<7,w^o) 
= - (8tt73<7 2 ) [l-g(0)], where g(0) is the pair- 
correlation function evaluated at zero distance; 

(iv) The following relations satisfied in the high- 
frequency limit by the real and imaginary parts of 
fxc for q<*>: \im a ^Ref xc (q,(o)=f xc (O,°°) + c/(0 3/2 
and lim^^ lmf xc (q,oo)= — c/co 312 , with c = 237r/15 
in the high-density limit evaluation of the irreduc- 
ible polarization propagator. 

These exact relations were used by Gross and Kohn 
(1985) to build a Pade approximation to the f xc (q 
^0,w). Although these results were obtained for the 
three-dimensional electron-gas case, similar relations 
were obtained for the two-dimensional case (Iwamoto, 
1984; Holas and Singwi, 1989) and extended to include 
two electron-hole pairs in both the longitudinal and the 
transverse response functions (Nifosf etal., 1998). 



Higher exact sum rules (up to seventh order) have been 
derived by Sturm (1995), and they can be used as a strin- 
gent test on the quality of the response functions derived 
from approximated f xc functional. 50 

Additional constraints on the matrix elements of the 
exchange-correlation kernel were obtained by Gonze 
and Scheffler (1999) using the Keldysh formalism 
adopted by van Leeuwen (1998). In particular, at the 
exact exchange level the following exact relation has to 
be satisfied close to the Kohn-Sham resonance energy 
wij : 

( < fyl/ XC (^/)l*«;) = <^l^-y x |^) 

-<^|s,-yj^)-<« ; -|yk>, 

(A5) 

where v x F is the Hartree-Fock nonlocal operator evalu- 
ated with the Kohn-Sham wave functions and V x is the 
exact exchange potential in DFT The last term describes 
an unscreened electron-hole attraction in contrast to the 
effective interaction in the TDLDA approximation (as 
discussed above). This is due to the fact that Eq. (A5) 
for the kernel explicitly shows two parts, namely, the 
first two contributions containing a shift linked to 
Hartree-Fock electron addition and removal energies, 
and the second one correcting with respect to those en- 
ergies, whereas the TDLDA kernel never passes explic- 
itly through electron addition and removal energies. In 
the single-pole approximation, expression (A5) leads to 
a correction of the excitation energies with respect to 
Kohn-Sham eigenvalue differences, which turns out to 
be identical to the correction obtained by Gorling (1996) 
in a first-order perturbation theory, in the difference be- 
tween the many-body and the second-quantized Kohn- 
Sham Hamiltonians. In addition, the inclusion of corre- 
lation effects leads to dynamically screened matrix 
elements through the same unscreened exchangelike 
term and an additional screened Coulomb interaction 
(Gonze and Scheffler, 1999). The similarities of this ap- 
proach to the GW quasiparticle correction and the in- 
clusion of electron-hole (excitonic) effects within the 
Bethe-Salpeter equation of many-body theory are de- 
scribed in Sec. VI. 

Another rigorous constraint is known as the 
harmonic-potential theorem (Kohn, 1961; Dobson, 
1994). It deals with the motion of an interacting many- 
electron system confined in a parabolic potential well 
under a spatially uniform time-dependent external field. 
This system exhibits sharp resonances at the bare 
harmonic-oscillator frequency, independently of the 
electron-electron interaction. This exact result for the 
response of an interacting electron system stems from 
the invariance of the harmonic potential under a trans- 
formation to a homogeneously accelerated reference 
frame (Vignale, 1995). In this case the dynamics of the 
electronic center of mass is completely decoupled from 



50 In particular, the TDLDA and RPA satisfy the first three 
odd-frequency sum rules but not the seventh (Sturm, 1995). 
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that of the internal degrees of freedom. This puts inter- 
esting constraints on approximate theories of time- 
dependent many-body physics such as TDDFT. Note 
that the TDLDA approximation to TDDFT satisfies the 
harmonic-potential theorem because the exchange- 
correlation potential follows the density when the latter 
is moved (that is, the exchange-correlation potential is 
local in both time and space; see below). The Gross and 
Kohn (1985) approximation violates the harmonic- 
potential theorem as shown by Dobson (1994), but it is 
possible to perform a simple modification for frequency- 
dependent local theories to satisfy this theorem (Vignale 
and Kohn, 1996; Dobson, Biinner, and Gross, 1997; Vig- 
nale et al, 1997). The essential idea is to make the action 
functional depend on the relative density p re/ (r,f)=p(r 
+ R CM (t),t), where RcmU) is the time evolution of the 
center of mass of the system. In this spirit an approxi- 
mate frequency-dependent exchange-correlation func- 
tional for small displacements from the equilibrium den- 
sity p can be constructed, as was done earlier by Gross 
and Kohn (1985), but now satisfying the harmonic- 
potential theorem. Thus the exchange-correlation po- 
tential reads V xc {*,t) = V^ A {r)+ P t dt' f xc [p Q {r),t 

— t']Sp re i(r,t'), with f xc as given in Gross and Kohn 
(1985). 51 

From the previous general relations for f xc one can 
see that a local-density approximation for time- 
dependent linear response in general does not exist, as 
long as one keeps on describing dynamical exchange- 
correlation effects in terms of the density only. This 
means that f xc is a strongly nonlocal functional of the 
density. This problem can be overcome if the current 
density is added as a basic variable in a generalized 
Kohn-Sham scheme (Vignale and Kohn, 1996). In this 
way one can derive a local current-density-functional 
theory of both current and density responses valid for 
slowly varying densities and potentials. This scheme was 
generalized by Vignale et al. (1997) to include dynamical 
contributions beyond the linear response. In particular, 
exchange and correlation beyond the TDLDA lead to 
the appearance of complex and frequency-dependent 
viscoelasticity stresses that, in the homogeneous 
electron-gas case, provide an additional damping mecha- 
nism to the decay into electron-hole pairs. 

The fact that the kernel is very nonlocal is also con- 
firmed by the calculation of the correlation energy of the 



A different perspective in functional development is 
achieved by the fact that the virial theorem for the exchange- 
correlation potential has been shown to hold for time- 
dependent electronic systems (Hessler et al., 1999). Moreover, 
the time dependence of the exchange-correlation energy is 
solely determined by the exchange-correlation potential 
[dE xc /dt =jdr dp(t,t)l dt V AC (rf)]. These exact relations have 
important implications for the construction of approximate 
time-dependent functional. In particular, the TDLDA ap- 
proximation is shown to fail badly in regions where the time- 
dependent density differs considerably from its ground-state 
counterpart (Hessler et al, 1999). 



homogeneous electron gas by Lein et al. (2000). They 
have shown that the nonzero spatial range of f xc {r,r',a>) 
cannot be neglected, whereas the frequency dependence 
is less important as far as total correlation energies are 
concerned. In fact, Gonze et al. (1997) have indicated 
that/ vc must have a 1/q 2 divergence, for g— >0. This con- 
tribution is supposed to solve the "metal-insulator para- 
dox," i.e., the seeming contradiction concerning systems 
that are in reality nonmetallic, but have a vanishing 
Kohn-Sham gap. 

APPENDIX B: DERIVATION OF THE EQUATIONS 
COMMON TO THE TDDFT AND BETHE-SALPETER 
APPROACHES 

1 . Dyson-like equation for the macroscopic dielectric 
function 



Given a matrix of the form 



M~- 



m 



III) 



m , 



m 2 m 

with m 00 being a c number, its inverse is 
\ 1 



(Bl) 



m 



-i 



+ 



X 



{m,QQ—m\m l m 2 ) 

T -1 

— m l m 



m 1 m 2 m l m 2 m\m 1 



(B2) 



If M is the dielectric function e = l — vP in its matrix 
form in reciprocal space (G, G'), then 



' M~ 



1 

'00 



: [ e oo e \ e le i\ 



The quantities e 00 , e\ , e 2 , and e are given by 



^OO" 



1 00' 



[s 2 ] G ' = -i>g'-Pg'o> G'#0, 
[elGG'^GG'-^G^GG' , G,G'#0. 
This yields 

e M = 1 "^o^oo - 2 ^o-Pog^gg^g'-Pg'o- 



G,G'#0 



Now one defines 



GG' 



5 G.G' 



VgP 



G r GG' > 



(B3) 

(B4) 
(B5) 
(B6) 
(B7) 

(B8) 

(B9) 



with 



(B10) 



0, G=0 
v G , G^0' 

i.e., v indicates that the G = contribution is to be left 
out in the bare Coulomb interaction. In matrix notation 
this yields 

S-GG'^ ° T e ). (BH) 
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Then the inverse of e is given by using Eq. (B2), 
1 r \ 



GG' 



(B12) 



Thus Eq. (B8) leads to Eq. (2.9), if in the latter P is 
defined as 



E -PCG'+ 2 PGKB KK ,V K r(q)P K r G r 



^GG' — ^GG'' ^-i 1 U«"n' 
K,K' 

This gives the wanted matrix expression, 
P = P+Pe~ 1 vP, 



(B13) 



(B14) 



and, together with e 1 = (1 — vP) , the Dyson-like 
equation (2.10). 

2. Effective two-particle equations 

In order to solve a Dyson-like equation such as Eq. 
(2.20), one has to invert a four-point function for each 
frequency. This problem can be reformulated as an ef- 
fective eigenvalue problem; in fact, the physical picture 
of interacting electron-hole pairs suggests using a basis 
of eigenfunctions ip n of the effective one-particle Hamil- 
tonian, from which the starting density matrix had been 
constructed, expecting that only a limited number of 
electron-hole pairs will contribute to each excitation 
(see Fetter and Walecka, 1971). These functions are sup- 
posed to form an orthonormal and complete set. Any 
four-point function 5 can then be transformed as 

= 2 ^*j(ri) ^n 2 0l) ^3(12) ^ 4 ( r 2) 5 («i« 2 )(« 3 « 4 ) ■ 



(B15) 

The four-point independent-quasiparticle polarization 

"lQP> 

P '/e/ 3 ( r l> r 2' r 3> r 4) 

_ y, (fn-fn>)<fin( r l)'Pn>(r2)<t'*>(U)<Pn(r3) 

e n-en'-<* 

(B16) 

is then diagonal in this basis and 

reads PlQP (« 1 « 2 )(»3"4) = (fn 2 ~fn l )^n 1 ,n 3 ^n 1 ,n 4 I ( € n 2 

-e n -w). With the definition n = [l - 4 P 1QP Ky\ Eq. 
(2.20) becomes S = H 4 P]q P . In transition space, 

n( ni „ 2 )(„ 3 „ 4 ) = [(e m2 -e mi -a>)S mitm3 S m2!m4 

(/fflj /m 2 )^(m 1 m 2 )(m3m 4 )](« 1 « 2 )(« 3 /i 4 ) 

X(e„ 4 -e„ 3 -w). 
Defining an effective two-particle Hamiltonian, 

+ (/n 1 _ /« 2 )^(/i 1 « 2 )(n 3 « 4 ) > (B17) 



one can then rewrite S as 

5(„ 1 n 2 )(„3„ 4 ) = [^ 2p -/«](„ 1 1 „ 2 )(„ 3 „ 4 )(/„ 4 -/„ 3 )- (B18) 

The one-particle transition energies on the diagonal are 
the eigenvalues of the starting one-particle Hamiltonian. 
Schematically, the effective Hamiltonian [Eq. (B17)] has 
the form 



H 2p 

(«1»2)("3"4) 



A B 
D 



(B19) 



with 



A = 



(n±n 2 )i 






{v'c'} 


{CV} 


{vc} 
{cv} 


jj2p,res 
tl (vc)(v'c') 

~lK( vc )( c i v >)]* 


K( vc )( c > v >) 

rTj2p,res -, ^ 
\- n (vc)(v'c')i 




We'} 


(B20) 

wn 


B = {vc} 
{cv} 




K( vc )( c >z>) , 

~^(cv)(c'c') 



(B21) 



D= {vv} 
{cc} 



{v'v'} 



{c'c'} 



(^-f«)<5««"5««' 

' e;- e :: \<\ : /<\,> 



(B22) 



The resonant part is defined as 

H 2 ( P v 'c)( v ' c ') = { € c~ € v)8v,v'<>c,c' + K( uc )( v > c >- ) . (B23) 

It corresponds to transitions at positive absorption fre- 
quencies w. Due to the factor (/„ — /„ ) in Eq. (B18), 
only the first column of 

W (n 1 1 » 2 )(« 3 n 4 ) = [^ 2P - /w ](« 1 1 , !2 )(n 3 » 4 ) . ( B24 ) 

i.e., with (n 3 n 4 ) = {c'v '} and {v'c'}, contributes to the 
calculation of S. 

Defining A a = A — Io), it is clear that only A ~ 1 is going 
to be relevant, in other words, the part A in Eq. (B20). 
This means that only pairs containing one filled and one 
empty Bloch state contribute, which is physically mean- 
ingful. It reflects the fact that only the hole-particle part 
[Eq. (3.7)] of the two-particle Green's function yields the 
absorption and loss spectra, as pointed out in Sec. III. A. 

It is this part A in Eq. (B20) which is referred to as 
the two-particle Hamiltonian H 2p in the present work. It 
is in general not Hermitian and can be further separated 
into four blocks: two blocks on the diagonal with the 
transition energies and the interaction kernel K, and 
two off-diagonal coupling blocks with contributions only 
from the interaction kernel: 
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H 2p-_ 



j_^2p,re$ ^-coupling 
_ |~ ^coupling j * ^jj2p,res^ J ' 

The resonant part is Hermitian, ^fjZpses^* 



(B25) 



: ^fjZpses^ wn jj e coupling part alone is symmetric, 

^coupling = (^coupling) r Jhe ^ j n ^ , ower right j g de _ 

noted antiresonant. Neglecting the coupling part is 
called the Tamm-Dancoff approximation. 

The spectral representation of the inverse two-particle 
Hamiltonian is 



[H 2p I^(n in2 )(n 3 n^)-^J 
A , A 



A: 



\,k' \ f 



E x -ai 

(B26) 

which holds for a system of eigenvectors and eigenval- 
ues of a general (not necessarily Hermitian) matrix de- 
fined by H 2p A x = E x A x , where N KtK * 
= 2,, „ A*^" 1 " 2 ^!^ 1 is the overlap matrix of the gen- 
erally nonorthogonal eigenstates of H 2p . Hence H 2p is 
diagonalized, and from its eigenvalues E x and eigen- 
states A" 1 " 2 the four- and two-point quantities of inter- 
est are constructed. In particular, one obtains the mac- 
roscopic dielectric function given in Eq. (2.23), and an 
equivalent expression for x- 



APPENDIX C: AN f xc FROM THE BETHE-SALPETER 
APPROACH 

In this subsection, we derive a time-dependent 
density-functional formalism for excitation spectra from 
the many-body Bethe-Salpeter equation (Reining et al. , 
2002). The results obtained using an approximation to 
this kernel, denoted as RORO, were already discussed 
in Sec. V.C. 

The main idea is to compare the effective two-particle 
Hamiltonians H 2p - TDDFT and H 2p < BSE derived for TD- 
DFT and the Bethe-Salpeter equation, respectively. If 
the quasiparticle and Kohn-Sham eigenfunctions are 
equal, and if all matrix elements of h 2 p< tddft anc j 
H 2p ' BSE involving those transitions which actually con- 
tribute to the spectrum are equal, than the resulting 
spectra will be identical. This yields the condition 



/«(q,G,G')= 2 



1 



\n x ,n 2 ;G) 



n l n 2 n 3 n 4 (fn 1 fn 2 ) 

x^(„ 1 „ 2 ) ( , l3 , l4 )(1>*)" 1 («3,"4;G'). (B27) 

The factor (/„ — /„ ) can never be zero for the particle- 
hole and hole-particle contributions of a nonmetal. The 
matrices $ are defined as 



(B28) 



and the operator T is 



_QP_ ,QP_ -DFT I DFT 



" " -I- e \ r) /*) 

e - 2 «i ' °n 1 n i °n 2 n 4 



• r (n 1 n 2 )(n 3 n 4 ) \ t n, c n 1 



(B29) 



It is now clear that, if the transformation from real to 
transition space x^ijj n (xi) was complete in all four in- 
dices, Eq. (B27) could never be satisfied, since otherwise 

H 2p,TDDFT an( j H 2p.BSE should alsQ be equal j n real 

space — and they cannot be, due to the way the S func- 
tions are put. On the other hand, if the two operators 
cannot be made equal, then the spectra can be equal 
only if at least one of the two operators (in that case,/^) 
is energy dependent. However, in practice only a finite 
number of transitions contribute to the optical spectrum. 
This means that one can use an incomplete basis in tran- 
sition space. Therefore one can still find a static operator 
that satisfies the required equality in transition space in 
a particular energy range, even though the real-space 
operators are not equal. (The invertibility of the matri- 
ces <t> may, however, be questionable in some particular 
cases). 

In view of the ongoing discussions about the 
exchange-correlation kernel it is interesting to examine 
some of its features, and in particular its long-range be- 
havior for solids. To this end, we note that 
<t> v c ,k+q,k(G = 0) goes to zero as q for small q. Since 
W vc vc behaves as a constant, if <t> is invertible this im- 
plies immediately that / xc (q,G=G' =0) behaves as 
1/q 2 . Note that there are two such long-range terms 
coming from (a) the electron-hole attraction (of nega- 
tive sign) and (b) the energy shift between the quasipar- 
ticle and the DFT eigenvalues [of positive sign, as pre- 
dicted by Gonze et al, 1997, on the basis of their study 
of the polarization-dependence of the exchange- 
correlation energy (Gonze et al, 1995; Resta, 1996), and 
as verified by a comparison of Kohn-Sham and experi- 
mental linear and nonlinear susceptibilities by Aulbur 
et al, 1996]. A discussion of how this approximated ker- 
nel reproduces the optical spectrum of semiconductors 
for the case of bulk silicon is presented in Sec. V.C. 
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